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signals  are  generated  using  both  the  HARPO  and  UMPE  models.  The  arrival 
structures  are  investigated  and  the  relation  between  features  in  the  arrival  structures 
for  matching  and  the  physical  parameters  are  identified.  Some  insight  into  the 
performance  of  the  SIFT  algorithm  is  gained  which  relates  matching  and  physical 
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I.  INTRODUCTION 

Acoustic  source  localization  has  been  an  intensive  research  issue  for  the  past  few 
decades.  Prior  to  that,  passive  SONAR  was  used  to  obtain  an  estimate  of  the  source  direction 
and  other  techniques,  such  as  Target  Motion  Analysis  (TMA),  Contact  Motion  Analysis 
(CMA)  and  triangulation,  were  used  to  get  an  estimate  of  the  range  to  the  source.  With  the 
introduction  of  faster  computers,  other  possibilities  arose.  One  of  the  techniques  that  was 
developed  is  known  as  Matched  Field  Processing  (MFP).  The  MFP  process  consists  of 
systematically  placing  a  synthetic  source  at  each  point  in  a  search  grid  and  comparing  the 
signal  received  from  all  the  synthetic  sources  with  the  signal  received  from  the  true  source. 
When  the  synthetic  source  location  and  the  true  source  location  match,  the  correlation  of  the 
true  and  the  synthetic  received  signals  should  be  maximum.  Most  of  the  work  in  this  field 
has  been  done  with  array  receivers  and  narrow  band  signals.  The  project-in-hand  concerns 
itself  with  the  case  of  a  single  receiver  hydrophone  and  broadband  signals. 

The  objective  of  this  study  is  an  examination  of  the  influence  of  the  physics 
mismatch  in  the  prediction  of  the  acoustic  propagation  on  a  number  of  MFP  algorithms  for 
a  single  hydrophone  receiver  and  a  transient-like  point  source.  A  time-domain  signal 
autocorrelation  matching  is  performed  to  produce  the  ambiguity  surface  for  the  localization. 
A  full  wave,  parabolic  equation  model  is  employed  to  produce  a  synthetic  signal  and  a 
reciprocal  prediction  to  provide  a  baseline  for  the  correlation.  Predictions  from  a  ray-based 
propagation  model  are  then  matched  to  the  synthetic  signal.  By  comparing  this  ambiguity 
surface  with  the  baseline,  the  influence  of  the  physics  mismatch  due  to  the  ray  predictions 
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can  be  quantified.  Several  aspects  of  the  ray  model  can  be  independently  affected  providing 
information  of  model  degradation  on  localization  performance.  A  secondary  but  not  less 
important  objective  is  to  suggest  directions  for  future  research  beyond  the  present  available 
algorithms. 

The  remainder  of  this  thesis  consists  of  five  chapters.  Chapter  II  describes  the 
propagation  models  and  arrival  structure  synthesis.  Chapter  III  describes  the  matching 
algorithms  and  the  influence  of  physical  parameters  on  these  algorithms.  Chapter  IV 
describes  the  setup  of  the  experiments  and  the  results  of  the  autocorrelation  matching  and 
the  approximate  autocorrelation  matching.  Chapter  V  describes  the  issues  which  were 
encountered  during  the  research  and  were  not  investigated,  but  which  may  lead  to  more 
insight  into  the  performance  of  present  algorithms  or  lead  to  more  robust  algorithms. 
Chapter  VI  presents  the  conclusions  of  the  study.  In  the  Appendix,  a  derivation  for  a 
frequency  based  analogue  of  a  generalized  beamformer  for  a  single  hydrophone  is  presented, 
and  the  relationship  to  autocorrelation  matching  is  briefly  discussed. 


H.  PROPAGATION  MODELS  AND  ARRIVAL  STRUCTURE  SYNTHESIS 

Propagation  models  and  in  particular  the  relation  of  their  parameters  to  physical 
parameters  form  a  major  part  of  the  basis  of  this  research.  The  step  from  the  output  of  the 
propagation  models  to  an  arrival  structure  from  another  important  part.  This  is  also  the  step 
where  some  of  the  major  approximations  are  made.  For  this  research,  two  propagation 
models  have  been  used:  a  full  wave  parabolic  equation  model  and  a  Hamiltonian  raytracing 
propagation  model.  The  models  are  based  on  the  same  wave  equation  but  different 
assumptions  are  made  to  generate  an  approximate  solution.  We  begin  with  the  homogeneous 
wave  equation  for  the  acoustic  pressure  in  a  medium  with  sound  speed  c{x)  and  density 
p(J0,  (Jensen  etal.  1994) 


p(x)V 


I   p(*) 
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Assuming  a  harmonic  solution,  p(x,t)=p{5t)e™1 ',  leads  to  the  corresponding  Helmholtz 
equation 
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where  S  is  the  source  function  and/?  is  pressure,  which  are  generally  complex.  The  right- 
hand  term  is  zero  when  the  source  is  not  included,  which  is  true  for  the  major  part  of  the 


ocean  For  a  fixed  density  and  a  range  independent  environment  the  formulas  could  even  be 
more  simplified  to 
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p(f)=0.  (3) 


Although  the  models  incorporate  more  features  and  are  based  on  generalized  versions  of  the 
above  formulas  (e.g.,  including  currents)  those  features  are  not  used  in  this  project. 

A.        HAMILTONIAN  RAYTRACING 

Raytracing  or  geometrical  acoustics  assumes  that  sound  waves  travel  along 
geometric  paths  called  rays.  Each  'ray'  travels  along  a  certain  path  which  is  dictated  by  the 
local  sound  speed  and,  where  appropriate,  boundary  interactions.  Assuming  a  constant 
density,  we  assume  a  solution  of  the  form  (Jensen  et  al.,  1994) 

.    ,    -    A  (x) 

/=o  (/o>y 

which  represents  a  ray  series  solution.  This  form,  although  in  general  divergent,  provides  an 
asymptotic  approximation  to  the  exact  solution  of  the  Helmholtz  equation  in  certain  cases. 
The  series  term  expresses  the  frequency  dependence  of  the  amplitude. 


Usually  only  the  first  term  in  the  series  solution  is  used  in  the  derivation.  The 
governing  equations  for  the  ray  theory  are  then 
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and 

2VT(^).V4.(f)+V2T(fM,  =  -V24(^)     7  =  1,2, ,  (7) 

where  ris  the  travel  time  phase  and  A  is  the  amplitude.  The  first  equation  is  known  as  the 
eikonal  equation.  It  determines  the  local  direction  of  the  phase  front  and,  as  rays  travel 
perpendicular  to  phase  fronts,  it  determines  the  ray  paths.  The  second  and  latter  equations 
are  called  the  transport  equations  of  which  only  the  first  is  commonly  used.  This  is  the  same 
as  taking  a  first  order  approximation  instead  of  a  series  solution.  The  second  and  latter 
equations  determine  the  amplitude  along  the  ray  path.  Furthermore,  the  expressions  above 
make  use  of  the  approximations 
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such  that  variations  of  the  sound  speed  and  variations  in  amplitude  occur  over  much  larger 
scales  than  a  wavelength.  The  neglect  of  higher  order  terms  in  the  above  solution  leads  to 
the  assumption  that  the  energy  is  conserved  within  a  ray  tube  (higher  order  terms  include 
leakage)  and  leads  to  a  break  down  of  the  transport  equation  at  or  near  focal  points 
(caustics). 

In  geometric  acoustics,  a  ray  behaves  much  like  a  particle.  Therefore,  Hamiltonian 
theory  can  be  applied  to  the  propagation  of  sound  (Lighthill,  1978).  The  first  order  ray 
equations  derived  from  the  eikonal  equation  can  be  written  as 

dx 

—=k(xylk(x),  (9) 
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where  s  is  the  path  length,  k(x)  = is  the  wave  number,  k0= —  is  the  reference  wave 

c{x)  c0 

number,  c0  is  the  reference  sound  speed  and  n  is  the  acoustical  index  of  refraction.  The 
equivalent  Hamiltonian  equations  can  be  written  as 

dx 

—  =Vj:H(x,k,u,t),  (12) 
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ax 


H(x,k,u,t)=uz-c2(x)k2=0,  (14) 


and 

T  =  4^  <15> 

In  this  information,  the  Hamiltonian  is  zero  along  the  ray  path.  For  this  project,  time  and 
frequency  dependence  of  the  Hamiltonian  are  neglected  (i.e.,  no  volume  attenuation, 
currents,  or  changing  environments).  All  frequency  dependent  features  such  as  boundary 
interactions  are  calculated  in  the  postprocessing.  Although  the  Hamiltonian  is  a  constant  of 
the  motion  in  this  formulation,  the  number  of  degrees  of  freedom  (3)  exceeds  the  number 
of  constants.  In  general,  ray  paths  are  then  chaotic  (Smith  et  al.,  1992).  However,  in  this 
project,  the  environment  is  two  dimensional  and  range-independent,  and  all  the  ray  paths  are 
regular. 

B.         HAMILTONIAN  RAYTRACING  PROGRAM  FOR  THE  OCEAN  (HARPO) 

The  raytracing  program  used  for  this  thesis  is  the  Hamiltonian  Raytracing  Program 
for  the  Ocean  (HARPO).  HARPO  is  written  in  Fortran  and  its  roots  go  back  to  a  program 


written  by  Dudrak  (1961).  It  has  since  been  updated  by  many  different  people.  The  program 
treats  the  ocean  as  a  continuous  changing  environment  and  not  as  a  large  number  of 
stratified  segments  thereby  avoiding  effects  such  as  false  caustics.  Of  the  ray  parameters, 
only  the  ray  paths  and  the  travel  times  are  computed  by  the  HARPO  program.  Amplitude  and 
path  length  are  determined  in  the  postprocessing. 

The  heart  of  the  program  is  a  fourth-fifth  order  Adams  Moulton  predictor  corrector 
algorithm  to  integrate  the  Hamiltonian  equations  along  the  ray  path.  The  accuracy  of  the 
integration  process  and  the  spacing  of  the  initial  ray  fan  are  extremely  important  parameters 
for  the  calculation  of  the  eigenray  parameters.  Comparisons  with  analytical  results  of  travel 
times  have  indicated  that  the  calculations  are  highly  accurate  with  phase  errors  £  n/100  up 
to  10  km  at  1  kHz. 

Two  features  in  the  design  of  the  program  make  the  standard  output  not  particularly 
suitable  for  a  straightforward  eigenray  extraction  routine.  The  rays  are  referenced  to  a 
spherical  coordinate  system  and  the  points  along  the  ray  are  not  equally  spaced  in  range. 
This  hampers  some  of  the  necessary  improvements  for  eigenray  extractions  close  to  the 
boundaries,  and  makes  it  necessary  to  use  a  local  Cartesian  coordinate  system  to  avoid 
complicated  formulations.  In  what  follows,  a  two-dimensional  environment,  specifying 
range  and  depth  positions  (r,z)  will  be  assumed. 


C.        EIGENRAY  EXTRACTION 

Eigenrays  are  rays  that  connect  source  and  receiver.  The  number  of  eigenrays  is 
generally  quite  large.  Usually  only  a  small  number  has  enough  amplitude  to  be  important, 
however.  Eigenray  extraction  programs  are,  in  a  general  sense,  root  solvers.  In  Fig.  1,  a  plot 
of  launch  angle  versus  depth  at  some  finite  range,  it  is  clearly  shown  that  the  eigenrays  for 
a  given  range  and  depth  are  the  roots  for  that  depth  of  the  curve  given  in  the  figure. 


Figure  1    Launch  angle  versus  depth  plot  at 
a  certain  range. 


As  the  output  of  HARPO  is  not  organized  in  a  convenient  range/depth  grid,  the  process  of 
finding  the  eigenrays  is  more  involved.  An  eigenray  extraction  program  developed  in  Matlab 
(ray3d.m)  by  Chiu  deals  with  this  matter.  The  program  selects  the  rays  with  a  (local) 
minimum  vertical  distance  to  the  receiver  location  out  of  the  total  fan  of  rays  and  calculates 


the  ray  parameters  by  a  third  order  polynomial  interpolation  using  these  selected  rays  and 
adjacent  rays.  No  retracing  of  the  eigenrays  is  done.  The  eigenray  parameters  calculated  by 
the  program  are:  phase  travel  time  r,  the  path  length  s,  the  local  sound  speed  c,  the  local 
sound  speed  gradient  Vc,  the  local  direction  of  the  ray  6,  the  launch  angle  60,  the  amplitude 
due  to  reflections,  the  phase  shift  due  to  reflections  and  turning  points,  and  the  amplitude 
due  to  geometrical  spreading. 

The  amplitude  along  the  ray  is  not  calculated  in  HARPO,  but  in  the  eigenray 
extraction  routine  based  on  the  equivalent  WKB  formulation  of  the  ray  solution.  The 
calculation  is  based  on  the  conservation  of  power  within  a  ray  tube.  The  formula  based  on 
the  assumption  of  a  point  source,  used  for  the  amplitude  calculation  is  (Ziomek,  1995) 

,                 .    oAr^dr^z)             /^cos(6„) 
Pirjf^ir.j.f    Ko -2 KISL ,  (16) 

where  \Pft-^)\  is  the  pressure  amplitude  at  the  gridpoint,  \PJ^r^x)\  is  the  pressure  amplitude 
at  1  m  (which  can  be  related  to  the  source  level),  p/r,z)  is  the  density  at  the  gridpoint, 
Po(ri>zi)  is  the  density  at  1  m  from  the  source,  c(r,z)  is  the  sound  speed  at  the  gridpoint, 
c(rltZj)  is  the  reference  sound  speed  at  1  m,  R0=\  m,  fi0  is  the  launch  angle  of  the  eigenray, 
r  is  the  horizontal  range  from  the  source  of  the  gridpoint,  and  dr±  relates  to  the  distance 
perpendicular  to  the  direction  of  the  wavefront.  We  can  rewrite  this  using 
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Po(ri^i)=Po(r'z)'  (17) 


*0  =  1> 


6=rsin((3), 


and  arrive  at 


_  c{r^z)     cos(P0) 
a'dr^)    r6/d$Q'  (18) 


This  formulation  is  known  to  produce  infinite  amplitudes  at  and  near  turning  and 
focal  points.  However,  the  numerical  implementation  used  here  prevents  the  amplitude  from 
blowing  up.  Ordinarily  the  d/dp0\s  calculated  using  the  selected  ray  (which  is  locally  closest 
to  the  gridpoint).  When  this  ray  is  too  close,  the  fan  is  opened  up  a  little.  The  same  is  done 
when  the  amplitude  exceeds  the  cylindrical  spreading  amplitude. 

The  routines  for  the  reflection  coefficients  are  based  on  plane  wave  reflection 
coefficients.  The  surface  reflection  coefficient  is  set  to  unity,  with  a  180  degrees  phase  shift, 
Rs=e~'*,  while  the  bottom  reflection  coefficient  is  calculated  including  the  sediment 
attenuation  a  (dB/Hz/km)  by 
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c  .=— a^-20log10(e), 
c''    2ti/    1000  10 


(19) 


0,=cos 


V  b    b>1  cos(6.) 


(20) 


so 


V 


P^cfe+^,/)sin(e/)"Pcsin(6z>) 
Pb(cb+icb  ,)sin(0,)+pcsin(0fc) 


(21) 


Written  in  terms  of  an  amplitude  and  a  phase  this  leads  to 


^=wA 


(22) 


where 


l*J  = 


(p^sin(0t)  -pcsin(db))2  +(pbcb  ,sin(0  .))2 
\|  (P^6sin(0.)+pcsin(e6))2+(p6^.sin(e.))2 


(23) 


and 


4>,=tan 


Pbcb,M%) 


-tan 


p,c    sin(0;) 


p6cfrsin(0.)  -pcsin(6fc)  P^sin(0 .) +pcsin(0fc) 


(24) 
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At  ray  turning  points,  we  can  express  the  phase  change  in  terms  of  an  equivalent  "reflection 
coefficient",  i.e.,R=e  ~tn/2  .  The  total  phase  for  selected  rays  (which  are  representative  of 
the  eigenrays)  is  then  calculated  as 

e  =  -tf,*/2-tf,*-X>  (25) 

where  the  number  of  turning  points,  surface  reflections  and  bottom  reflections  are  given  by 
N„  Ns  and  Nb  respectively. 

The  phase  of  a  ray  changes  according  to  reflections,  refractions  and  travel  time. 
Refractions  close  to  the  surface  and  close  to  the  bottom  may  not  cause  a  discrete  -te/2  phase 
change  but  something  between  -n/2  and  the  phase  change  due  to  the  reflection.  Furthermore, 
the  position  in  range  is  not  strictly  localized  at  the  refracted  points.  This  together  with  the 
fact  that  small  errors  accumulate  per  bounce  and  may  give  the  phase  of  a  ray  a  random 
component  which  can  be  quite  large.     . 

The  estimation  of  the  amplitude  of  a  ray  is  also  numerically  difficult.  Opening  up  the 
ray  fan  as  discussed  before  does  not  give  you  a  true  estimate  but  rather  a  conservative 
estimate  due  to  the  averaging  caused  by  opening  up  the  ray  fan.  However,  travel  times  are 
estimated  very  accurately.  Still,  at  high  frequencies,  small  errors  may  lead  to  unacceptable 
phase  errors.  For  ranges  and  frequencies  considered  in  this  project,  we  do  not  expect  travel 
time  errors  to  be  significant. 
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D.         ARRIVAL  STRUCTURE  SYNTHESIS  FROM  EIGENRAY  DATA 

Several  different  approaches  have  been  used  to  derive  arrival  structures  from 
eigenray  data  for  different  localization  algorithms.  The  important  parameters  used  in  the 
synthesis  are 

Tn  eigenray  travel  time  (relative  and  absolute  travel  time), 

6  eigenray  phase  (due  to  reflections  and  refractions), 

fc  carrier  frequency  (which  together  with  the  travel  time  may  account 

for  a  major  part  of  the  phase), 
a„  eigenray  amplitude,  and 

s(t),S(f)  source  function  (pulse/transient)  in  time  or  frequency  domain  (real 

or  complex). 
For  the  synthesis,  the  ocean  is  considered  as  a  linear  time  invariant  filter  causing  a  time 
delay  and  a  phase  shift.  The  total  amplitude  and  phase  contribution  due  to  the  reflections  and 
refractions  for  each  eigenray  can  be  written  as 


W- 


-V, 


IL*jn 


n*jf) 


ri  tun 


7=1 


(26) 
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where  the  product  in  Eq.  (25)  represents  the  accumulated  effects  of  multiple  reflections  and 
refractions  and  the  last  equation  specifies  the  conjugation  for  negative  frequencies.  For  this 


9(/)={ 


project  we  assume  that  Rn  is  frequency  independent.  The  frequency  response  of  a  real 
received  signal  from  an  individual  arrival  can  then  be  written  as 

YJfrHJft  X{f)  (28) 

where 

Hn{f)=an{f)  Rn  Tn(f)  (29) 

is  the  complex  frequency  domain  transfer  function,  a  „  represents  the  geometrical  amplitude 
factor,  and  T„  represents  the  phase  change  due  to  the  time  delay.  Note  that  Yn(J)=Y*(-f) 
since  the  time  domain  signal  is  real.  We  assume  that  the  amplitude  factor  is  also  frequency 
independent  and  the  delay  factor  can  be  written  as 

Tn(f)=e-i2*K  (30) 
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We  can  now  write  YJf)  as 


YJf)=an  Rn  e^'e-^'Xtf)-  (31) 


In  the  time  domain,  this  can  be  written  as 


-id     i2nAt-i„) 


yjf)°°.  k.  \x^    'e       J#  02) 


The  arrival  structure  is  then  defined  as  the  sum  over  all  eigenrays  of  the  individual  arrivals 


i.e. 


N 


y®=EyJto-  (33) 


n  =  \ 


For  some  applications,  it  may  be  easier  to  first  sum  in  the  frequency  domain, 

*</>=£  YJf),  (34) 


and  then  transform  the  results  to  the  time  domain, 


M=fYV)ea«df.  (35) 


These  two  approaches  should  be  equivalent. 
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For  some  localization  algorithms,  it  may  be  preferable  to  perform  all  the  calculations 
in  the  frequency  domain.  For  bandpass  signals,  it  might  be  desirable  to  cast  it  in  terms  of  the 
preenveiope  or  complex  envelope  of  the  received  signal.  This  results  in  a  few  minor  changes 
in  the  formulas.  Specifically,  we  treat  the  time  domain  signal  as  complex  (i.e.,  neglect  the 
negative  portion  of  the  spectrum)  and  define  the  basebanded  frequency  response  as 

Yn(f)=Hn(f+fc)  Xif)  (36) 

or 

f,(fi=a,  Rn  M  e  ■'%  -2nf-z-e  -W".  (37) 

If  we  write  the  complex  source  function  as 

x(t)=F-x{X(f)}  (38) 

where 

x{t)  =  {x{t)+ix{t))e-iW<>  (39) 

then  the  complex  received  signal  of  a  single  eigenray  can  be  written  as 


y(t)=x(t-Tn)e 

=f(r-Tn)cos(e+27i/;Tw)-/f(/-Tn)sin(e+27c/cT„) 

where/,  is  the  carrier  frequency  or  center  frequency  of  the  pass  band. 


(40) 
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Another  way  to  write  the  received  signal,  which  was  used  in  previous  simulations 
and  which  is  especially  applicable  to  low  pass  signals,  can  be  derived  starting  from  the 
preenvelope  (no  base  banding)  of  the  source  function 

xp(t)=x(t)+ix(t\  (41) 

where  x(t)  is  the  Hilbert  transform  of  the  real  source  function  x(t)y  and  the  frequency 
spectrum  of  the  preenvelope  x  (t)  is  the  same  as  that  of  x(t)  but  with  the  negative  frequencies 
discarded.  From  the  preenvelope  of  the  received  signal 

y^xfi-xje-*",  (42) 

the  real  received  signal  can  be  described  by 


yn{i)=Re{xp{t-xn)e~lQ" 

^(r-T„)cos(6„)+ix(r-Tw)sin(e„). 


(43) 


E.         THE  PARABOLIC  EQUATION  MODEL 

The  parabolic  equation  model  used  here  is  based  on  an  approximation  to  the 
Helmholtz  description  of  the  wave  equation  in  a  cylindrical  coordinate  system.  The  majority 
of  the  ocean  environment  is  well  suited  for  a  description  in  cylindrical  coordinates. 
Assuming  a  time  harmonic  solution,  the  Helmholtz  equation  in  cylindrical  coordinates  takes 
the  form  (see,  for  example,  Jensen  et  al.,  1994) 
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I  ± 

r   Br 


'   dpp-j$)\  i    &pfr,z,$)    cPpfrxQ)     2  2/      xx   ,      as     „    d*,-  -x     ,^ 

r— i- — J— + — J— +k0n2(r,z,<$>)p£r,z,<b)=-4nP0b(r-rs),    (44) 

V         or       )  rz         d$z 


dz- 


where 


pirj,§fy)=pfrj:,$)e 


2%ft 


(45) 


G) 


kQ= —  is  a  reference  wave  number  and  the  acoustic  index  of  refraction  is  defined  by 


n=- 


c(rj:,§) 


(46) 


As  energy  is  primarily  radiating  outward  from  the  source,  p/r,z,(f>)  can  be 
approximated  by 


p/rf,$)=V/rxWfaorl 


(47) 


where  ¥/r,z,  <j>)  is  a  slowly  varying  envelope  and  H0  (k0r)  is  the  outward  going  Hankel 
function.  In  the  far  field,  the  Hankel  function  takes  the  form  of  an  outward  going  plane  wave 
and  so  we  define 


p/r,z,(t))=-lT(^(j))e"'V 


(48) 
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Substituting  this  into  Eq.  (43)  yields 


2m        ~>2 


d2x¥  ...  aT    i   a2Y    a 

+/2A:n + + + 

dr2  dr     r2    d<|>2    dz2 


■  2r„2 


*b(«    "1  + 


4Lr: 


Y=0, 


(49) 


where  the  source  function  on  the  right-hand  side  has  been  dropped  Neglecting  the  azimuthal 
coupling  and  the  near  field  terms  and  assuming  that  Fis  a  slowly  varying  function  with 
range,  we  arrive  at 


Defining  the  operators 


and 


this  can  be  written  as 


ai=JL  #1^  (bM)<F. 

dr     2A0    Sz2      2 


(SO) 


op      2 


~lkQ     


(51) 


op  «\\  /» 


(52) 


'Vf^v^- 


(53) 
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With  the  operators  defined  above,  this  constitutes  what  is  commonly  referred  to  as  the 
"standard"  parabolic  equation  (Tappert,  1977).  For  this  work,  we  have  employed  the  higher 
order  "wide  angle"  parabolic  equation  (Thomson  and  Chapman,  1983)  with  operators 
defined  by 


rp         _  dr 

1 WAPE~' 


dz 


_2 


1+ 1     +1 

dz: 


(54) 


and 

UWAPE=-(n-l),  (55) 

The  parabolic  equation  model  that  was  used  for  this  project  is  the  University  of 
Miami  Parabolic  Equation  model  (Smith  and  Tappert,  1994).  A  split-step  Fourier  algorithm 
is  used  to  numerically  integrate  the  solution  in  range.  This  involves  alternatively  applying 
the  Uop  and  the  Top  operators  in  the  z-domain  and  the  1^-domain,  respectively,  where  each 
operator  is  simply  a  scalar  multiplier.  The  algorithm  for  stepping  in  range  from  r  to  r  +  Ar 
can  then  be  succintly  expressed  as 

T(r+A^)=e'^^(^xFFr{^-'^f-Wx[FFr(T*(^))]*},  (56) 
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where  the  wide  angle  f    operator  in  the  kz-space  is  defined  as 


op 


TwAPE^J-}- 


1- 


1/2 

(57) 


The  output  of  the  model  is  in  the  form  of  the  field  functions  ^(magnitude  and  phase)  and 
has  been  referenced  to  a  source  strength  of  1  uP  at  1  m.  The  field  is  only  computed  at 
discrete  points  and  the  spacing  in  depth  and  range  are  the  primary  parameters  that  determine 
the  accuracy  of  the  results.  The  wavelength  of  interest  in  this  study  is  1.5-2  m  and  the  grid 
spacing  used  is  40  cm  in  depth  and  50  cm  in  range.  Note  that  the  grid  size  is  ~  X/4  which  is 
necessary  to  obtain  the  highest  accuracy  up  to  ±  90  degrees  propagation  angles.  In  general, 
we  would  expect  that  a  courser  grid  size  would  produce  similar  results  in  shallow  water 
environments  where  propagation  angles  are  often  limited  to  <  15  degrees. 

F.         PE  ARRIVAL  STRUCTURES 

The  pressure  is  defined  in  terms  of  the  PE  field  function  by 

p(r*A=-±R#***(r*A  (58) 

ft 

where  P0  is  the  amplitude  of  the  source  pressure  measured  at  R0=l  m.  The  broadband  results 
were  obtained  by  running  the  model  multiple  times  for  all  discrete  frequencies  in  the 
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bandwidth  under  consideration  and  then  performing  a  Fourier  synthesis  to  yield  travel  time 
results.  Assuming  a  normalized  source  amplitude,  the  complex  arrival  structure  of  the 
pressure  field  can  then  be  written  as 

p(r^t)=fp(r^J)ei2^dt 

—  oo 

(59> 

=-Lfeik<rx¥(z,rJ)ei2nf'df. 
yfrL 


ilnf— 

Since  e '  ^  -e      c° ,  this  phase  factor  can  be  neglected  by  defining  the  reduced  time  T-t- — 
such  that 


p{r^T) =-1  (¥  (zs,4>)e  i2^df.  (60) 


Note  that  using  a  reduced  time  has  no  influence  on  the  autocorrelation  function. 
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HI.  LOCALIZATION  ALGORITHMS 

The  Transient  Localization  Project  at  the  Naval  Postgraduate  School  has  in  recent 
years  studied  different  localization  algorithms  for  the  scenario  of  one  receiver  hydrophone 
and  a  point  source.  All  the  routines  were  based  on  fundamental  concepts  of  generalized 
correlation  functions  described  in  Bendat  and  Piersol,  1996  and  most  of  them  are  described 
in  Miller  et  al.,  1996. 

Localization  algorithms  may  be  considered  generalized  beamformers  in  which  the 
plane  wave  replicas  have  been  replaced  by  more  complicated  replicas  of  the  acoustic 
propagation  (e.g.,  modes,  beams,  or  the  vertical  pressure  field).  The  algorithms,  usually 
referred  to  as  processors,  are  in  most  cases  based  on  a  Hermitian  quadratic  product.  The 
exact  form  is  determined  by  the  constraints  that  are  put  on  the  processor  output. 

The  main  algorithm  developed  previously  developed  at  the  Naval  Postgraduate 
School  was  the  Signal  Integration  Filtering  in  Time  (SIFT)  algorithm  which  is  a  form  of  time 
domain  autocorrelation  matching  (Benson,  1995).  This  algorithm  had  considerable  success 
in  the  localization  experiments  in  the  Barents  Sea  (Benson,  1995).  Although  the  results  were 
favorable,  there  were  still  questions  about  the  influence  of  environmental  model  parameters 
and  acoustic  (ray)  model  parameters  on  the  results. 
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All  the  localization  algorithms  described  below  are  described  for  the  single 
hydrophone  case  but  can  easily  be  extended  to  an  array  of  hydrophones.  The  source  is 
described  as  an  omnidirectional  point  source  and  the  receiver  is  a  single  omnidirectional 
hydrophone. 

All  localization  algorithms  are  based  on  matching  the  received  signal  with  a  replica. 
Replicas  are  simulated  received  signals  of  synthetic  sources  at  a  large  number  of  gridpoints 
in  a  search  space.  When  the  propagation  model  is  100%  correct  and  the  synthetic  and  real 
source  functions  are  equal,  the  real  and  synthetic  received  signal  should  exactly  match  when 
the  real  and  synthetic  source  positions  coincide.  Results  of  the  localization  are  usually 
represented  as  an  ambiguity  surface.  All  MFP  algorithms  are  based  on  this  principle  in  some 
form.  Using  the  reciprocity  principle  reduces  the  amount  of  work  to  a  manageable  size,  and 
makes  localization  possible  in  a  reasonable  and  operational  feasible  time.  The  reciprocity 
principle  states  that  in  an  environment  without  time  variations  (e.g.,  currents)  the  acoustic 
pressure  at  location  B  from  an  omnidirectional  source  at  location  A  and  the  acoustic 
pressure  at  A  from  an  equivalent  source  at  B  are  inversely  proportional  to  the  densities  at 
A  and  B.  If  we  assume  the  densities  are  the  same  at  A  and  B,  we  may  interchange 
hydrophone  and  source  locations  and  thereby  reduce  the  amount  of  work. 
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A.        FULLY  COHERENT  LOCALIZATION 

Fully  coherent  localization  is  based  on  a  normalized  time  domain  cross  correlation 
of  a  signal  and  the  set  of  replicas.  Given  a  received  signal 


r(t)=fR{f)ei2^df, 


(61) 


and  predicted  replicas 


p{t)=fP(f)ei2^df, 


(62) 


the  cross  correlation  is  given  by 


fR*(f)P(f)ei2^df 


CM 


\ 


CO  oo 


f\R(f)\2dff\P(J)\2df 

— oo  — oo 


(63) 


In  the  frequency  domain,  the  replicas  may  be  represented  by 


P(f)=H{f)S(f) 


(64) 
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where  H(f)  represents  the  propagation  model  and  S(f)  represents  the  source  function  of  the 
replica,  or  template  signal. 

The  maximum  of  the  peak  will  reflect  the  degree  in  which  the  modeled  signal 
matches  the  real  received  signal.  Only  the  peak  value  of  this  cross  correlation  is  utilized. 
When  the  synthetic  and  real  source  coincide  and  the  ocean  model  is  100%  correct  the  value 
should  be  1.  This  method  suffers  from  a  number  of  flaws.  R(f)  and  P(f)  are  in  general 
complex  and  their  phase  is  very  dependent  on  the  ocean  model  parameters.  (We  neglect  the 
source  function  dependence  which  will  be  discussed  in  the  chapter  on  future  work).  A  phase 
mismatch  of  90  degrees  at  a  given  point  will  produce  a  zero  cross-correlation.  Such  a 
mismatch  at  1000  Hz  can  easily  be  caused  by  small  sound  speed  or  bathymetry  errors,  so  this 
method  does  not  seem  to  be  viable  for  practical  situations.  Applying  the  algorithm  at 
basebanded  signals  does  not  solve  the  phase  problem. 

B.         SEMI-COHERENT  LOCALIZATION 

To  avoid  the  problems  with  phase  mismatch,  which  are  more  likely  at  high 
frequencies,  a  semi  coherent  approach  can  be  taken.  In  this  approach,  the  absolute  value  of 
the  time  domain  replicas  and  the  received  signal  is  passed  through  a  lowpass  filter.  This  is 
similar  to  broadband  demodulation.  The  matching  algorithm  is  then  performed  on  the  low 
pass  signal,  i.e., 
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R'fcHJfiflfWye'2*'^  -»*dt  (65) 


P'fcHuMflfPV^eW'dfle  -i2^dt  (66) 


-oo   ~oo 


and  so 


^R'(f)P'if)efl%fxdf 


\ 


(67) 


f\RW4ff\PW4f 


The  match  is  again  found  at  the  maximum  of  C^r).  Previously  this  algorithm  was  used 
without  much  success  (Miller  et  al.  1996)  and  will  not  be  used  here.  However,  careful 
analysis  of  the  effects  of  errors  in  the  parameters  of  the  environment  or  the  eigenrays  on  the 
arrival  structures  may  shed  some  light  on  the  poor  results. 

The  replica/* '(/)  shows  that  the  largest  amplitudes  coincide  with  the  earliest  arrivals. 
This  portion  of  the  signal  therefore  provides  the  largest  contribution  to  the  result  of  the 
match.  However  the  amplitudes  of  these  arrivals  are  most  prone  to  error  due  to  poor 
estimation  of  the  amplitudes  of  refracted  rays  which  form  the  largest  part  of  the  initial 
arrivals.  The  semi-coherent  algorithm  may  avoid  the  phase  errors  but  amplifies  the 
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amplitude  errors.  An  improvement  could  be  to  normalize  the  amplitude  of  the  replica.  A 
number  of  other  improvements  could  also  be  made,  such  as  sampling  at  the  arrival  times  of 
the  individual  rays,  or  clipping  the  template. 


C.         TIME-DOMAIN  AUTOCORRELATION  MATCHING 

One  of  the  problems  with  the  previous  algorithms  is  the  lack  of  an  absolute  time 
reference.  The  time  domain  autocorrelation  matching  algorithm  removes  the  absolute  time 
reference  from  the  received  signal  and  the  replica.  The  autocorrelation  of  the  received  signal 
is  defined  as 


GJix^fR'WWe^df  (68) 


and  the  autocorrelation  of  the  replica  is  similarly 


Gpp(i)= JP  *(f)P(f)e  i2^df.  (69) 
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The  autocorrelation  matching  can  then  be  computed  by 


E  GRR(x)Gpp(z) 


Arp~ 


\ 


E  i<^)i2E  \GPP^t 


(70) 


or  in  the  frequency  domain,  using 


G^fG^e-^df, 


(71) 


GPP(f)-jGpp(x)e-'2^df 


(72) 


then 


E  Gmif)Gppif) 


^RP~' 


\ 


E  \Gn0T,  \GpP(f)r 


(73) 


Originally  the  matching  was  designed  for  the  real  signals  in  the  time  domain. 

There  are  two  primary  advantages  of  autocorrelation  matching.  First  the  influence 
of  noise  tends  to  concentrate  near  zero  lag  values.  By  cutting  out  values  around  zero  lag,  the 
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effect  of  noise  is  significantly  reduced.  Cutting  out  values  around  zero  lag  also  has  a  distinct 
effect  on  the  contrast  of  the  ambiguity  surface  as  can  be  seen  in  Fig.  2.  Matching  values  now 
range  from  -1  to  1.  Values  between  -1  and  0  are  neglected  in  all  subsequent  plots.  The 
second  advantage  is  the  removal  of  any  absolute  time  of  reference.  Only  relative  arrival 
times  of  separate  multipaths  are  compared. 

Throughout  most  of  this  project,  the  received  signals  will  be  generated  synthetically 
using  the  PE  model.  Replicas  were  generated  using  the  same  PE  model  for  various 
environments  to  determine  the  optimum  number  of  lag  values  near  zero  to  remove.  The 
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Figure  2   Effect  of  notching  out  near  zero  lag  values  of  autocorrelation  matching  at  a 
certain  range.  Plots  (a)  -  (f)  show  increasing  number  of  notched  out  values. 
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results  indicate  that  this  number  depends  on  the  shape  of  the  source  signal  and  therefore  on 
the  bandwidth.  Furthermore,  where  it  has  a  positive  effect  on  the  contrast,  it  has  a  negative 
effect  on  the  footprint  size.  For  example,  by  removing  seven  points  near  zero  lag  (zero  lag 
plus  three  at  positive  and  three  at  negative  lags)  for  a  negative  gradient  environment,  the 
footprint  is  reduced  from  about  five  meters  in  depth  to  about  two  and  a  half  meters  in  depth. 
The  footprint  size  does  not  seem  to  be  reduced  further  when  more  points  than  the  optimum 
number  of  points  are  removed.  However,  the  removal  of  more  points  begins  to  degrade  the 
performance  by  increasing  the  level  of  false  peaks  (sidelobes).  Note  that  the  removal  of  near 
zero  lag  points  in  the  time  domain  is  equivalent  to  removing  the  mean  and  lowest  order 
trends  from  the  frequency  domain.  We  found  that  the  removal  of  points  corresponding  to  a 
width  of  about  2/bandwidth  is  adequate  under  most  circumstances  for  source  signals  having 
a  Blackman  shaped  frequency  response. 

D.         SEMI-COHERENT  AUTOCORRELATION  MATCHING 

As  has  been  done  for  semi-coherent  localization  a  less  coherent  form  of  localization 
can  be  defined  for  the  autocorrelation  matching.  Using  Eq.  71  and  Eq.  72  we  can  define  a 
semi-coherent  autocorrelation  matching  algorithm  as 
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This  form  reduces  the  influence  of  the  phase.  Working  with  only  positive  values,  however, 
makes  the  dynamic  range  of  the  ambiguity  surface  very  small.  Results  for  PE  generated 
synthetic  received  signals  and  PE  generated  replicas  look  promising  but  the  performance 
when  comparing  PE  received  signals  and  ray  model  replicas  was  not  adequate  to  pursue  this 
algorithm  at  this  time. 

E.         THE  SIFT  LOCALIZATION  ALGORITHM 

The  SIFT  (Signal  Integration  Filtering  in  Time)  algorithm  has  been  developed  as  a 
result  of  previous  research  (Miller  et  al.  1996,  Benson  1995).  It  is  based  on  a  ray  code 
solution.  The  templates  are  generated  using  the  arrival  time  and  phase  of  each  eigenray. 
These  parameters  determine  uniquely  each  gridpoint  in  a  search  grid.  The  SIFT  algorithm 
assumes  that  the  signal  can  be  written  as  the  sum  of  scaled,  phase  and  time  shifted  versions 
of  the  original  signal,  i.e., 
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where  a„  is  the  real  amplitude  and  xn  and  0n  are  the  travel  time  and  phase,  respectively,  of 
an  eigenray,  x(t)  is  the  real  source  signal  and  x{i)\s  its  Hilbert  transform.  The  autocorrelation 
of  the  replicas  is  given  by 
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The  correlation  between  individual  eigenrays  can  be  written  as 


or 


G  (x)=aacosQcosQR(x+x-x) 
+a  a  cos6  sin6  /?U(t+t„  -t  ) 
+a  a  sin0  cos©/?- (t+t  -x) 
+a  a  sin0  sin0  R.Xx+x  -x  ). 
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Using  properties  of  the  Hilbert  transform,  this  can  be  rewritten 

G     (x)=a  a  G  (t+t  -t  )cos(0  -0  ) 
+a  a  G  (t+t  -t  )sin(0  -0  ), 

n    m     xx^-  n       m'         v     m        «-" 


(79) 


where  G    is  the  Hilbert  transform  of  the  real  source  autocorrelation,  G^ 

XX  ** 

In  the  original  SIFT  algorithm,  three  major  approximations  are  made: 

1)  the  amplitudes  are  set  equal  to  unity  (a=l  for  all  n); 

2)  the  phases  are  independent  of  frequency; 

3)  the  autocorrelation  functions,   G^  andG^  ,  are  approximated  by  discrete 
functions  defined  by 


G  (t+t  -t  )=<         T  "T"  T"  (80) 
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and 
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0      otherwise 


where  6t  is  chosen  to  be  the  correlation  time  of  the  measured  signal. 
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The  matching  part  of  the  algorithm  is  then  defined  by 
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where  Rjj  and  R^  are  the  vectors  of  the  sampled  correlations  of  the  replicas  and  the  received 
signal 
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Usually  the  autocorrelation  of  the  received  signals  will  be  formed  as  a  linear  correlation  and 
will  then  be  sampled.  The  above  description  assumes  the  eigenrays  are  sorted  in  order  of 
increasing  travel  time. 

A  similar  derivation  can  be  made  for  the  representation  of  a  complex  signal.  The  ray 
arrival  structure  can  then  be  represented  by 

■Z"Jllh>^  (85) 

w  =  l 

N 


»=1 


In  general,  we  may  want  to  use  signals  that  have  been  basebanded  or  shifted  in  frequency 
by  some  amount^.  The  complex  source  may  then  be  defined  by  replacing 
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x(t-z)=*x(t-x)e  ~""*B\  (86) 


The  autocorrelation  of  the  received  signal  can  now  be  expressed  as 


Gf=fy(t)y*(t+z) 
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To  simplify  the  expression  we  define 
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and  then  the  correlation  between  the  signals  of  the  individual  eigenrays  can  be  expanded  as 
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or 
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If  we  would  apply  similar  approximations  as  with  the  original  SIFT  algorithm  we  would 
approximate  the  amplitude  with  unity  and  the  autocorrelation  function  would  have  to  be 
sampled  with  a  single  point  approximation  for  T=xm-Tn.  When  taking  only  the  positive  or 
only  the  negative  lags  into  account  for  the  approximation,  the  matching  function  would  be 
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When  both  the  negative  and  positive  lags  are  used  the  result  will  obviously  be  real.  For 
bandpass  signals  this  seems  a  worthwhile  approach.  An  alternative  approach  to 
accommodate  bandpass  signals  would  be  to  shift  them  to  the  lowest  possible  frequency 
band,  form  the  real  autocorrelation  of  the  received  signal  and  then  apply  the  original  SIFT 
algorithm  for  real  signals. 

F.         SIFT  AND  THE  PHYSICS  OF  THE  LOCALIZATION  PROBLEM 

The  choice  of  St  is  one  of  the  major  difficulties  with  the  SIFT  algorithm.  One  way 
to  optimize  the  algorithm  would  be  to  try  and  optimize  the  results  for  a  small  range  of  6t's. 
For  simulations  6t  can  easily  be  calculated  from  the  cross  correlation  of  x(t)  and  x(t).  The 
performance  may  be  enhanced  by  optimizing  over  a  small  range  of  6t's,  although  this 
reduces  the  computational  speed.  For  this  approach  to  be  valid,  it  is  also  necessary  that  dt 
is  less  than  the  separation  between  distinct  arrivals  and  greater  than  the  time  sampling  rate. 
For  very  broadband  signals  with  good  time  resolution,  this  approximation  may  be  very  good. 
However,  this  method  is  quite  crude  and  may  not  work  for  narrow  band  signals.  An 
alternative  may  be  to  use  a  multi-point  approximation  instead  of  a  three-point 
approximation. 

The  3-point  approximation  of  the  correlation  is  a  low  frequency  approximation.  At 
high  frequencies  and  for  bandpass  signals  this  approximation  is  no  longer  valid.  Furthermore 
the  values  with  which  the  approximation  is  made  (1,-1,1)  are  dependent  on  the  shape  of  the 
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source  function.  It  may  be  worthwhile  to  look  at  a  less  rigid  approximation.  A  few  test  cases 
were  examined  in  which  the  Hilbert  transform  of  the  autocorrelation  was  neglected  and  no 
serious  degradation  in  performance  was  observed. 

The  phase  is  one  of  the  most  sensitive  parameters  used  in  the  algorithm.  Small 
changes  in  phase  at  turning  points  and  boundary  interactions  give  the  phase  at  larger 
distances  a  more  or  less  random  value.  Although  arrival  structures  may  look  similar,  the 
autocorrelation  will  be  very  different,  and  much  of  the  match  is  lost.  For  the  complex 
envelope  version  of  SIFT  the  result  not  only  depends  on  the  phase  difference  due  to 
reflections  and  refractions  but  also  on  the  phase  difference  due  to  the  time  delay.  Although 
the  travel  times  are  known  to  be  a  very  robust  physical  parameter  in  these  cases,  small 
errors  may  become  significant  at  high  frequencies. 

In  general,  the  calculation  of  eigenray  amplitudes  is  reasonably  accurate.  However 
there  are  a  few  arrivals  in  those  environments  where  the  estimation  of  amplitudes  is  quite 
poor.  Typically,  these  rays  correspond  to  refracting  rays  near  caustics  producing  erroneously 
large  amplitudes.  A  large  arrival  will  dominate  the  autocorrelation  and  make  it  look  very 
much  like  the  arrival  structure.  This  may  be  a  good  reason  to  normalize  the  amplitude.  When 
you  normalize  the  amplitude,  the  later  arrivals  have  the  same  weight  in  the  correlation  and, 
therefore,  in  the  matching  algorithm  as  the  earlier  arrivals.  This  has  a  number  of  advantages. 
The  later  arrivals  carry  as  much  information  as  the  first  few  and,  therefore,  taking  them  into 
account  increases  the  amount  of  information  in  your  match.  Also  solutions  of  the  eigenray 
extraction  routine  are  more  stable  for  the  later  arrivals  than  for  earlier  arrivals. 
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IV.  EXPERIMENTAL  SETUP  AND  RESULTS 

To  gain  insight  into  the  performance  of  the  autocorrelation  matching  algorithms  and 
approximate  autocorrelation  matching  algorithms,  several  computer  experiments  have  been 
done.  The  experiments  have  primarily  been  performed  for  a  frequency  band  of  750-1000  Hz. 
The  baseline  for  the  results  was  defined  as  the  results  of  the  autocorrelation  matching  of  a 
PE  generated  source  and  PE  generated  templates.  For  most  of  the  experiments  a  source 
range  of  5500  m,  a  source  depth  of  59.8  m  and  a  receiver  depth  of  40.2  m  have  been 
assumed.  For  the  environments  defined  below,  at  this  range  the  first  arrivals  arrive  at 
approximately  3.7  seconds  after  the  transmission  and  the  significant  part  of  the  arrival 
structure  is  about  0.4  seconds  in  length.  As  the  project  did  not  focus  on  the  frequency 
characteristics  of  the  source  function,  the  same  shape  has  been  assumed  for  the  source  signal 
and  the  templates.  Throughout  the  project  analysis  the  absolute  scale  of  the  pressure  field 
calculation  has  been  neglected  as  the  similarity  of  the  arrival  structures  is  emphasized  not 
the  absolute  values  of  the  pressure  field. 

A.        ENVIRONMENTS 

For  the  synthetic  experiments,  three  shallow  water,  range  independent  environments 
have  been  defined.  The  first  has  a  negative  sound  speed  gradient  (1500-1475  m/s),  the 
second  has  a  constant  sound  speed  of  1500  m/s,  while  the  third  has  a  positive  sound  speed 
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gradient  (1500-1501.5  m/s).  All  other  acoustic  parameters  are  fixed  for  the  three 
environments.  The  water  depth  is  100  m  and  the  bottom  parameters,  the  attenuation, 
density  and  sound  speed,  were  chosen  as  a=0.48  dB/(km  Hz),  pbottom=1900  kg/m3, 
0^^=1650  m/s,  respectively,  consistent  with  a  sand  bottom.  The  bottom  is  modeled  as  an 
infinite  half  space  in  both  propagation  models.  The  propagation  loss  for  the  environments 
at  a  frequency  of  875  Hz  is  shown  in  Fig.  3  (a),  (b),  and  (c).  These  plots  suggest  a  very  rapid 
variation  of  the  acoustic  pressure  with  respect  to  depth.  The  vertical  size  of  localization 
footprints  is  therefore  anticipated  to  be  small. 

The  third,  positive  gradient,  environment  has  only  been  used  to  generate  a  limited 
set  of  ACM  baseline  results.  It  was,  perhaps  wrongly,  anticipated  that  the  results  might  not 
significantly  differ  from  the  results  generated  for  the  iso-speed  environment,  and  therefore 
the  computational  effort  would  not  outweigh  the  additional  results.  While  this  appeared  to 
be  the  case  for  the  third  environment  (with  a  weak  positive  gradient),  the  results  for  the  first 
environment  (strong  negative  gradient)  were  significantly  different.  The  analysis  suggested 
that  most  of  the  algorithms  performed  well  for  the  second  environment  (iso-speed)  but  were 
at  their  performance  limits  for  the  first  environment.  It  is  therefore  likely  that  a  less  severe 
sound  speed  gradient  environment  could  have  given  more  insight  into  the  exact  performance 
limits  of  the  algorithms. 
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Figure  3    Transmission  loss  at  875  Hz  for  (a)  a  negative 
sound  speed  gradient  (1500-1475  m/s)  environment,  (b)  an  iso- 
speed  (1500  m/s)  environment,  and  (c)  a  positive  sound  speed 
gradient  (1500-1501.5  m/s)  environment. 
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B.         HAMBLTONIAN  RAYTRACING  PROGRAM  FOR  THE  OCEAN 

The  Hamiltonian  Raytracing  Program  for  the  Ocean  (HARPO)  has  a  number  of 
parameters  which  directly  influence  the  accuracy  of  the  results.  To  be  able  to  use  the  results 
for  this  application,  a  maximum  ray  launch  angle  separation  of  0.05  degrees  has  been 
determined.  Larger  spread  caused  unacceptable  inaccuracies  in  the  eigenray  extraction 
program.  The  program  has  been  set  to  output  all  the  points,  and  the  integration  accuracy  has 
been  set  to  W6.  Less  accurate  settings  would  degrade  the  performance  of  the  eigenray 
extraction  program.  Travel  times  and  eigenray  calculations  have  been  compared  with  results 
from  an  analytical  model  and  proved  to  be  very  accurate.  To  reduce  the  computational 
burden,  the  assumption  has  been  made  that  at  the  ranges  of  interest  the  attenuation  due  to 
bottom  reflections  would  be  high  enough  such  that  source  angles  larger  than  critical  could 
be  neglected.  This  reduces  the  number  of  eigenrays  to  about  40-50  per  grid  point  for  the 
negative  gradient  environment  for  a  range  of  about  5  km.  For  shorter  ranges,  source  angle 
limits  of  ±40  degrees  have  been  used  resulting  in  more  than  1 8  eigenrays  per  grid  point  for 
ranges  above  1  km.  This  set  up  resulted  in  a  ratio  of  the  smallest  to  the  largest  eigenray 
amplitude  per  grid  point  of  smaller  than  0.02  for  95%  of  the  environment.  The  grid  point 
spacing  of  the  eigenray  extraction  routine  has  been  set  to  5  m  in  range  and  1.4  m  in  depth. 
The  depth  spacing  does  not  meet  the  initial  minimum  requirement  of  five  grid  points  per 
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dimension  of  the  footprint  (12.5  m  x  1.2  m  at  5.5  km  for  a  source  depth  of  59.8  m  and  a 
receiver  depth  of  40.2  m).  It  is  a  compromise  due  to  the  computational  speed  of  the 
eigenrays  extraction  routine. 

During  the  experiments,  four  eigenray  extraction  problem  areas  were  indicated.  Very 
abrupt  change,  or  little  change,  together  with  noise  due  to  round  off  and  interpolation  errors 
and  using  too  large  angle  separations  creates  either  too  few  or  too  many  eigenrays  or 
'spurious'  eigenrays  at  these  points.  For  the  iso-velocity  environment,  a  separate  program  has 
been  used  that  is  based  on  an  analytical  solution.  This  allowed  a  little  more  flexibility  and 
considerably  reduced  the  run  time. 

C.         UNIVERSITY  OF  MIAMI  PARABOLIC  EQUATION  MODEL 

The  University  of  Miami  Parabolic  Equation  Model  (UMPE)  was  used  with  the 
wide-angle  PE  approximation,  as  described  in  Chap.  II.  The  maximum  computational  depth 
has  been  set  to  409.60  m.  A  depth  spacing  of  40  cm  leads  to  250  grid  points  within  the 
water  column,  and  4-5  grid  points  per  wavelength.  As  a  comparison,  only  70  grid  points  in 
depth  have  been  used  to  cover  the  water  column  for  the  ray  model.  The  grid  spacing  in  range 
has  been  set  to  50  cm  which  is  about  the  minimum  for  an  accurate  prediction  in  the 
environments  described  above.  Results  have  been  stored  in  2.5  m  and  5  m  range  increments 
for  the  negative  gradient  and  iso-speed/positive  gradient  environments,  repectively.  The 
frequency  is  sampled  from  750  Hz  to  1000  Hz  in  257  frequency  bins.  At  a  range  of  6500  m, 
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this  frequency  resolution  still  allows  the  calculation  of  the  correlation  as  the  Fourier 
transform  of  the  power  spectral  density  without  degradation  due  to  wrap-around. 

D.        RECIPROCITY 

In  the  experiments,  the  source  and  receiver  are  assumed  to  be  a  point  source  and  a 
point  receiver.  The  source  spectrum  is  shaped  as  a  Blackman  window.  This  avoids  large 
sidelobes  in  the  time  domain  arrival  structure  and  leads  to  a  pulse  length  of  approximately 
8  ms.  To  optimize  the  eigenray  extraction  routines,  a  large  number  of  reciprocity  checks  for 
different  source  and  receiver  locations  have  been  done.  The  autocorrelation  function  of  the 
arrival  structures  proved  to  be  an  excellent  tool  in  cases  where  the  arrival  structures  looked 
alike.  Small  differences  in  amplitude  and  especially  phase  led  to  large  differences  in  the 
autocorrelation  function.  This  suggested  that  localization  algorithms  based  on  the  time 
domain  autocorrelation  would  be  very  sensitive  to  these  parameters,  particularly  the  phase. 

The  PE-results  showed  very  good  reciprocity.  The  ray  code  results  showed  good 
reciprocity  except  at  certain  locations  where  erroneous  amplitude  calculations  degraded  part 
of  the  arrival  structure.  Looking  at  the  difference  between  the  ray  code  and  PE  arrival 
structures,  it  appears  that  the  largest  difference  can  be  found  in  the  first  few  arrivals.  The 
tails  of  the  arrivals  structures  are  very  similar.  This  suggests  that  the  ray  code  approximation 
for  refracting  rays  is  less  accurate  than  for  reflecting  rays.  Perturbation  analysis  of  the  phase 
and  amplitude  of  the  rays  suggests  that  the  arrival  structures  and  autocorrelation  are  more 


48 


susceptible  to  phase  errors  than  to  amplitude  errors.  Results  also  suggest  that  certain  parts 
of  the  autocorrelation  are  more  affected  than  others.  Very  short  lags  and  very  long  lags  seem 
to  suffer  more  than  intermediate  values.  Knowing  that  short  lag  values  are  notched  out  and 
large  lag  values  will  have  no  significant  contribution  to  the  match  suggests  that  the 
algorithms  should  at  least  have  some  degree  of  tolerance  toward  amplitude  and  phase  errors. 
This  suggestion  is  based  on  visual  analysis  of  a  large  number  of  results,  but  no  statistical 
analysis  has  been  performed. 

E.         MATCHING  ALGORITHMS 

Numerical  experiments  have  focused  on  autocorrelation  matching  and  the  SIFT 
matching  algorithm.  The  baseline  for  the  results  has  been  defined  by  the  PE  autocorrelation 
matching  over  the  full  250  Hz  bandwidth  using  a  Blackman  shaped  frequency  response.  The 
autocorrelation  matching  experiments  can  be  separated  into  their  variations  of  the  basic 
algorithm.  The  following  variations  have  been  specifically  examined:  1)  the  effect  of 
notching  out  points  near  zero  lag;  2)  the  effect  of  using  a  smaller  bandwidth,  which  was 
anticipated  to  give  a  larger  footprint  size;  3)  the  effect  of  splitting  up  the  frequency  band  into 
several  smaller  frequency  bands;  4)  the  effect  of  a  lower  frequency  resolution;  5)  the  effect 
of  raising  the  tail  of  the  autocorrelation  function,  i.e.  to  increase  the  weight  of  large  lag 
values  (corresponding  to  cross  terms  of  early  and  late  arrivals);  6)  the  effect  of  matching  the 
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absolute  values  of  the  autocorrelations  functions;  and  7)  the  effect  of  amplitude  and  phase 
perturbations  on  the  frequency  domain  arrival  structure. 

Specifically  for  the  ray  results,  additional  effects  have  been  examined:  1)  the  effect 
of  raising  the  tail  of  the  arrival  structure  to  increase  to  weight  of  the  later  arrivals;  2)  the 
effect  of  neglecting  the  first  few  arrivals,  as  the  amplitude  and  phase  of  these  arrivals  are 
suspected  to  be  the  cause  of  most  of  the  bad  results;  and  3)  the  effect  of  matching  the  arrival 
structures  directly,  in  both  complex  and  absolute  form. 

The  SIFT  algorithm  has  been  considered  in  four  different  ways:  1 )  a  form  where  the 
complex  signals  have  been  base-banded  and  a  1 -point  approximation  to  the  autocorrelation 
function  is  used;  2)  a  form  where  the  real  signals  are  shifted  in  frequency  to  a  center 
frequency  of  half  the  bandwidth  and  where  subsequently  a  modified  form  of  the  original 
SIFT  algorithm  is  employed;  3)  a  form  where  the  signals  are  shifted  and  then  a  zero  mean 
version  of  the  original  SIFT  algorithm  is  employed  which  results  in  a  multiple-point  match; 
4)  variations  on  the  original  low  frequency  implementation  of  the  SIFT  algorithm.  For  all 
the  variations,  notching  out  points  near  zero  lag,  neglecting  early  arrivals,  and  variations  of 
the  parameters  that  describe  the  approximation  to  the  correlation  functions  applied  in  SIFT 
have  been  considered. 

The  grid  spacing  of  the  templates  has  been  based  on  the  footprint  size  of  the  baseline 
results  and  computational  arguments.  As  the  main  area  of  interest  is  beyond  5  km,  the  results 
of  autocorrelation  matching  at  5.5  km  for  the  negative  gradient  environment  have  been  used 
for  determining  the  footprint.  The  footprint  size  is  approximately  1.2  m  in  depth  and  12.5 


50 


m  in  range.  With  a  norm  of  5  points  per  dimension  of  the  footprint  this  would  lead  to  a  grid 
spacing  of  20  cm  in  depth  and  2  m  in  range.  Analysis  at  other  ranges  suggests  the  footprint 
will  be  smaller  at  smaller  ranges  and  larger  at  larger  ranges.  The  excessive  run  time  of  the 
eigenray  extraction,  the  available  disc  storage  and  run  time  for  the  matching  algorithms  have 
lead  to  a  grid  spacing  of  2.5  m  in  range  and  40  cm  in  depth  for  PE  generated  templates  in 
the  negative  gradient  environment,  5  m  in  range  and  0.4  m  in  depth  for  PE  generated 
templates  in  the  iso-speed  and  positive  gradient  environments,  and  5  m  in  range  and  1.4  m 
in  depth  for  the  ray  code  generated  templates.  For  the  negative  gradient  environment,  the 
eigenray  results  have  only  been  calculated  for  1000-2000  m  and  4500-6500  m  due  to  the 
computational  load  of  the  MATLAB  based  eigenray  extraction  routines.  Attempts  to  apply 
the  recently  available  MATLAB  c-compiler  were  not  successful.  The  presently  available 
compiler  still  uses  many  standard  MATLAB  routines  which  prevent  generation  of  efficient 
code.  Also  the  nature  of  the  problem  which  relies  on  using  complex  numbers  degrades  the 
ability  to  generate  efficient  code. 

To  generate  synthetic  'received'  signals,  three  basic  algorithms  have  been  used:  an 
impulse  model  based  on  Eq.  43,  a  time  domain  model  based  on  Eq.  33,  and  a  frequency 
domain  model  based  on  Eq.  37.  Where  the  impulse  model  is  an  inherently  low  frequency 
model,  the  time  domain  model  is  suited  primarily  for  band  pass  signals.  A  small 
modification  to  the  original  time  domain  based  model  developed  by  Chiu  made  it  suitable 
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for  low  frequency  applications.  In  the  time  domain  model,  the  tails  of  the  pulse  are  removed 
which  is  a  major  advantage  in  testing  algorithms  where  overlap  of  pulses  is  a  crucial 
parameter. 

The  autocorrelation  of  the  synthetically  generated  received  signals  has  been 
computed  using  several  different  methods.  For  autocorrelation  matching  using  frequency 
based  data  as  generated  by  the  UMPE  model,  it  may  seem  that  the  transform  of  the  power 
spectral  density  may  be  the  most  practical  solution.  However  the  circular  convolution  does 
not  work  well  with  pulse  signals  where  the  pulse  length  is  more  than  half  the  transform  time 
base.  For  this  reason  the  autocorrelation  for  this  project  has  primarily  been  computed  from 
the  time  domain  arrival  structure.  Where  in  previous  work  an  'unbiased'  autocorrelation 
estimate  has  been  used,  for  this  work  the  *biased'  autocorrelation  estimate  has  been  used 
which  is  more  suited  for  pulse  shaped  signals.  As  the  time  base  is  determined  by  the 
frequency  resolution,  it  may  be  beneficial  in  future  experiments  to  artificially  stretch  the 
time  base  before  performing  the  autocorrelation.  Results  of  applying  phase  disturbances  to 
the  PE  results  suggested  that  keeping  the  time  domain  arrival  structure  centered  in  the  first 
half  of  the  time  base  can  not  always  be  accomplished.  While  theoretically  PE  phase  errors 
should  have  no  influence  on  the  autocorrelation  results,  there  will  be  a  degradation  in  the 
results  due  to  the  above.  Thus,  while  theoretically  the  time  domain  based  autocorrelation 
may  be  preferable,  in  practice  the  frequency  domain  based  calculations  are  more  robust. 

No  use  has  been  made  of  the  usual  performance  measures  such  as  peak  to  sidelobe 
level  or  mismatch  in  range  and  depth.  The  reason  is  that  these  measures  tend  to  put  a 
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number  to  a  match.  This  work  tries  to  explain  how  and  why  the  match  will  degrade.  It  has 
been  observed  that  a  match  tends  to  become  amplified  along  certain  lines  of  high 
correlation/energy.  Furthermore,  when  the  matching  parameters  exceed  certain  values, 
depending  on  the  type  of  matching,  environment,  range  and  bandwidth,  the  results  suggest 
the  single  match  will  explode  into  several  'matches'  spread  around  the  actual  source  location. 
This  still  gives  an  indication  of  the  source  location  but  it  will  not  give  a  single  point  match. 

F.         COHERENT  MATCHING  RESULTS 

Coherent  matching  has  been  performed  in  the  time  domain  and  in  the  frequency 
domain.  The  frequency  domain  approach  may  require  an  additional  parameter  to  allow  for 
time  shifts.  All  presented  analysis  uses  a  time  domain  approach.  Results  may  slightly  differ 
due  to  the  time  resolution  used.  Results  for  the  negative  gradient  environment  with  a  source 
at  5500  m  display  sidelobe  levels  comparable  to  autocorrelation  matching  with  an  optimum 
number  of  lag  values  notched  out.  The  footprint  size  for  the  same  setup  was  found  to  be 
larger  for  the  coherent  matching  than  for  the  autocorrelation  matching.  Coherent  matching 
results  using  ray  code  templates  and  PE  generated  source  signals  showed  some  phase  shift 
between  the  templates  and  the  source  signals.  In  general,  we  were  able  to  localize  correctly 
up  to  ranges  of  more  than  5  km.  Short  ranges  resulted  in  very  high  sidelobes,  especially  in 
the  iso-speed  environment.  No  significant  changes  have  been  made  to  the  basic  coherent 
matching  algorithm.  In  many  cases  where  the  autocorrelation  matching  algorithms  failed  to 
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localize,  the  coherent  matching  still  gave  reasonable  results.  Taking  into  account  a  phase 
shift  in  the  result  of  the  match  improved  the  performance  significantly.  There  are  however 
many  drawbacks  to  coherent  matching,  and  there  are  locations  where  we  found  that  many 
templates  and  the  received  signal  were  almost  alike  which  led  to  very  high  sidelobes. 
Discretized  phases  may  be  a  convenient  way  to  improve  the  coherent  matching  algorithm. 
In  addition  some  analysis  has  been  done  on  a  less  coherent  approach  using  the 
absolute  value  of  the  time  domain  arrival  structure.  The  dynamic  range  was  improved  by 
subtracting  the  average  over  a  certain  area  of  the  matching  resuls.  Altough  baseline  results 
showed  promising,  the  results  using  PE  generated  received  signals  and  ray  code  generated 
templates  showed  very  high  sidelobes  which  prohibited  correct  localization.  A  similar 
approach  will  be  discussed  further  in  Chapter  V. 
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Figure  4   Coherent  matching  for  the 
negative  gradient  environment 
(baseline  result),  with  a  source  at 
5500  m  range  and  59.8  m  depth  and  a 
receiver  at  40.2  m  depth. 


Figure  5   Coherent  matching  for  the 
negative  gradient  environment  using 
a  PE  generated  source  signal  and  ray 
code  templates  with  a  source  at  5500 
m  range  and  5  9.8  m  depth  and  a 
receiver  at  40.2  m  depth.  Accounting 
for  the  phase  shift  results  in 
correct  localization. 
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G.        AUTOCORRELATION  MATCHING  RESULTS  (BASELINE) 

The  autocorrelation  matching  algorithm  is  described  in  Chap,  m,  Section  C.  To 
improve  upon  the  results  of  the  basic  algorithm,  the  following  variations  have  been 
considered:  notching  out  near  zero  lag  values;  enhancing  the  dynamic  range  by  subtracting 
the  average  over  a  certain  area  of  the  match  and  renormalizing;  increasing  the  significance 
of  large  lage  value  by  filtering;  changing  the  sample,  rate,  bandwidth  and  frequency 
resolution;  and  splitting  the  frequency  band  into  several  small  bands.  To  get  an  initial 
assessment  of  the  robustness,  amplitude  and  phase  perturbations  have  been  added  to  the  PE 
generated  received  signals  and  PE  generated  templates. 

Results  for  both  the  negative  and  positive  gradient  environments  and  a  source  at 
5500  m  range  suggest  that  notching  out  zero  lag  and  the  first  2-4  lag  values  (corresponding 
to  (1  or  2)/BW)  gives  an  optimal  performance  in  terms  of  peak  to  sidelobe  level.  For  the 
isospeed  environment  the  optimum  was  found  by  notching  out  four  times  the  previous 
number  of  lag  values,  having  still  larger  sidelobes  than  the  other  two  environments. 
Presumably  arrival  structures  are  more  similar  for  the  iso-speed  environment.  Results 
suggest  that  the  optimum  number  of  notched  out  lag  values  depends  on  range  and,  to  a  lesser 
extent,  on  bandwidth  (except  for  low  frequencies  where  the  bandwidth  dominates).  The  first 
few  lag  samples  represent  the  low  frequency  trends  in  the  autocorrelation.  They  are  more 
related  to  the  power  in  the  signal  rather  than  the  shape  and  will  dominate  the  match  when 
they  are  not  removed  Notching  out  these  lag  values,  however,  decreases  the  footprint  size. 
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Notching  out  more  than  the  optimum  number  of  lags  does  not  significantly  reduce  the 
footprint  size  any  further  for  the  negative  gradient  environment.  For  the  iso-speed 
environment,  however,  the  footprint  size  continues  to  decrease  further  when  the  number  of 
notched  out  lag  values  increases. 

As  an  alternative  to  notching  out  values  around  zero  lag  in  the  time  domain,  two 
approaches  have  been  considered.  Removing  the  mean  from  the  frequency  domain  results 
leads  to  similar  results  as  notching  out  only  the  zero  lag  value.  Removing  lower  order  trends 
would  require  filtering  of  the  power  spectrum,  which  makes  it  much  more  complicated  than 
notching  out  values  in  the  time  domain,  with  theoretically  the  same  results.  The  second 
approach  aims  to  increase  the  dynamic  range,  and  is  based  on  the  subtraction  of  the 
matching  value,  averaged  over  an  area,  from  the  individual  matching  values.  Contrary  to 
notching  out  values  around  zero  lag,  this  approach  has  no  real  physical  interpretation.  The 
results  are  comparable.  However  sidelobe  levels  are  slightly  higher,  and  it  does  not  have  the 
degrees  of  freedom  to  optimize  the  results  as  with  notching  out  lag  values  in  the  time 
domain. 

The  footprint  was  found  to  depend  on  range,  receiver/source  depth,  sound  speed 
profile,  position  of  the  source  relative  to  the  'lines  of  high  correlation/energy1,  and  on 
bandwidth.  It  was  found  for  the  cases  considered  that  the  general  trend  is  for  the  footprint 
to  increase  with  range.  Footprints  for  the  same  source/receiver  depth  were  found  to  be  much 
larger  than  for  significantly  separated  source  receiver  depths.  In  general,  footprints  tended 
to  be  larger  for  the  environment  with  a  small  positive  gradient  and    the  iso-speed 
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environment  than  for  the  environment  with  a  large  negative  gradient.  The  position  relative 
to  a  line  of  high  correlation/energy  determines  if  and  how  the  footprint  will  stretch.  The 
most  significant  difference  between  the  slight  positive  gradient  or  iso-speed  environment  on 
one  hand  and  the  severe  negative  gradient  environment  on  the  other  is  the  vertical  extent  of 
the  footprint  which  is  nearly  twice  as  large  for  the  slight  positive  gradient  and  even  larger 
for  the  iso-speed  environment  than  for  the  negative  gradient  case.  The  smaller  footprints  of 
the  negative  gradient  environment  can  be  directly  related  to  the  more  severe  interference 
patterns.  The  significantly  larger  footprints  for  the  same  source  receiver  depth  are  thought 
to  be  caused  more  by  energy  matching  than  by  matching  of  the  arrival  structure  shape. 

The  number  of  points  or  sample  rate  did  not  appear  to  be  a  crucial  parameter.  There 
are  some  results  suggest  that  a  sample  rate  of  1000  Hz  might  give  better  results  and  a  slightly 
larger  footprint  than  a  sample  rate  of  250  Hz.  For  this  reason  the  bulk  of  the  analysis  has 
been  performed  for  a  sample  rate  of  250  Hz  for  the  negative  gradient  environment  and  1000 
Hz  sample  rate  for  the  iso-speed  environment.  Contrary  to  our  initial  expectations,  smaller 
bandwidths  generally  produced  smaller  footprints  and  higher  sidelobes  than  larger 
bandwidths.  This  may  be  due  to  fewer  arrivals  being  sampled  at  smaller  bandwidths  cause 
the  the  autocorrelation  (and  its  match)  to  drop  off  rapidly  away  from  the  source  location. 
Larger  bandwidths  on  the  other  hand,  resolve  more  multipath  arrivals  which  may  broaden 
the  autocorrelation  and  generate  a  larger  footprint. 
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A  lower  frequency  resolution  showed  a  larger  footprint  but  also  large  sidelobe 
levels.  Although  the  number  of  cases  examined  is  limited,  the  results  suggest  that  larger 
bandwidths  are  preferable,  as  expected.  Results  also  suggest  that  it  may  not  be  necessary  to 
apply  any  window.  Results  show  that  the  degradation  due  to  sidelobes/ringing  in  the  time 
domain  may  be  outweighed  by  the  increase  in  performance  due  to  the  larger  bandwidth. 

Separating  the  750-1000  Hz  band  into  4  smaller  bands  and  averaging  the  results  led 
to  a  slightly  larger  footprint  as  compared  to  the  full  bandwidth  results  but  with  a  significant 
increase  in  sidelobe  levels.  Later  arrivals  are  more  separated  in  time  than  earlier  arrivals  and 
are  also  more  stable  in  these  environments  as  they  are  not  related  to  refracted  propagation 
paths.  To  be  able  to  exploit  them,  the  later  arrivals  should  be  enhanced.  In  the  PE  model  the 
individual  arrivals  cannot  be  manipulated.  Therefore,  an  alternative  approach  was  used, 
enhancing  the  tail  of  the  autocorrelation.  Results  suggest  that  indeed  some  improvement  may 
be  gained  with  this  approach  in  terms  of  the  sidelobe  levels,  especially  if  it  is  applied  to  both 
the  received  signals  and  the  templates.  The  footprint,  however,  appears  to  be  smaller  than 
for  the  full  bandwidth  case  (7.5  m  in  range  and  2.5  m  in  diameter  at  5.5  km  for  the  negative 
gradient  environment).  The  results  were  not  found  to  be  significant  enough  to  pursue  this 
approach  any  further.  It  may  be  that  in  more  complicated  environments  this  approach  should 
be  reconsidered. 
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Figure  6   Baseline  autocorrelation  matching  results  for  the 
negative  gradient  environment.  The  source  is  located  at  5500  m 
range  and  59.8  m  depth  and  the  receiver  at  40.2  m  depth.  Three 
notched  out  near  zero  lag  values  have  been  used,  (a)  Full 
range,  (b)  Range  expanded  near  the  source  location. 
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Figure  7    Baseline  autocorrelation  matching  results  for  the 
iso-speed  environment.  The  source  is  located  at  5500  m  range 
and  59.8  m  depth  and  the  receiver  at  40.2  m  depth.  Three 
notched  out  near  zero  lag  values  have  been  used,  (a)  Full 
range,  (b)  Range  expanded  near  the  source  location. 
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Figure  8   Baseline  autocorrelation  matching  results  for  the 
positive  sound  speed  gradient  environment.  The  source  is 
located  at  5000  m  range  and  40.2  m  depth  and  the  receiver  at 
59.8  m  depth.  No  window  has  been  applied  and  3  near  zero  lag 
values  have  been  notched,  (a)  Full  range.   (b)  Range  expanded 
near  the  source  location. 
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To  get  an  assessment  of  the  robustness  of  the  algorithms  amplitude  and  phase 
perturbations  have  been  applied.  The  effect  of  phase  perturbations  has  already  been 
discussed.  The  amplitude  perturbation  has  been  applied  as  a  uniform  random  fluctuation. 
Using  a  random  distributed  amplitude  fluctuation  up  to  +/-  40  %  only  showed  a  reduction 
in  the  peak  to  sidelobe  level,  but  no  significant  increase  in  sidelobe  levels  or  change  in  the 
footprint  size  were  observed.  This  indicated  that  the  algorithm  is  fairly  robust  for  amplitude 
fluctuations. 

H.        AUTOCORRELATION  MATCHING  (NUMERICAL  RESULTS) 

When  matching  PE  source  signals  and  ray  code  templates  for  the  negative  gradient 
environment  using  a  windowed  frequency  response,  results  appear  to  be  changing 
periodically  for  ranges  beyond  about  5  km.  Moving  the  source  over  a  range  increment  of 
4500-6500  m,  a  correct  localization  was  obtained  about  50%  of  the  time.  Of  the  remaining 
half,  some  regions,  did  not  show  any  localization  while  in  other  places  the  localization  was 
obscured  by  high  sidelobes.  Short  ranges  give  more  consistent  results  but  also  contain 
inherently  high  sidelobes.  Even  with  a  severely  inadequate  ray  set,  localization  could  be 
performed  at  ranges  up  to  2000  m. 
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Several  changes  have  been  made  to  obtain  more  consistent  results  at  longer  ranges.  Notching 
out  a  large  number  of  lag  values  sometimes  improves  the  results  but  not  consistently. 
Applying  a  window  that  amplifies  the  larger  lag  values  of  the  correlation  gave  some 
improvement  when  applied  to  both  the  received  signal  and  the  templates. 

From  the  good  localizations,  a  1.2%  difference  in  range  between  the  ray  code 
solution  and  the  PE  solution  was  found  for  the  negative  gradient  environment.  This 
difference  is  probably  due  to  phase  errors  in  the  PE  approximation  or  a  sub-optimal  choice 
for  the  reference  sound  speed  and  using  below  optimal  values  for  the  range  and  depth  grid 
spacing.  The  footprint  size  is  usually  very  small,  5-15  m  in  range  and  1.4  m  (1  grid  point) 
in  depth.  Neglecting  the  initial  arrivals  in  the  ray  code  solution  did  not  improve  the  results. 
It  also  seems  that  higher  sampling  rates  give  slightly  higher  matching  peak  values  but  do  not 
significantly  improve  the  results.  Splitting  the  frequency  band  into  several  smaller  bands  can 
improve  the  results  for  a  certain  set  of  matching  parameters,  but  the  results  are  not 
predictable.  No  solid  reason  could  be  found  explaining  the  changing  behavior  with  range. 
It  is  speculated  that  this  is  not  directly  related  to  the  bandwidth  (smaller  bandwidths 
sometimes  produced  better  results  than  larger  bandwidths)  but  may  be  related  to  the  modal 
interference  pattern,  the  phase  errors  in  the  initial  arrivals,  and  possibly  under-sampling. 

No  correlation  was  found  between  eigenray  amplitude  problem  areas,  eigenray 
extraction  problem  areas,  arrival  times,  arrival  time  differences,  depth  of  the  perigee  of 
refracted  rays,  or  number  of  refracted  rays  and  the  changing  behavior  with  range.  The  phase 
of  the  'half  match  (using  only  positive  or  negative  lag  values)  was  observed  to  vary 
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periodically.  But  as  it  is  composed  of  many  parts,  it  is  more  indicative  than  a  real  measure, 
and  it  shows  no  sudden  changes  for  locations  where  a  match  was  expected  but  not  found. 
A  fixed  phase  difference  in  the  time  domain  should  not  show  in  the  autocorrelation.  Still, 
as  other  methods  do  not  seem  to  suffer  as  much,  the  reason  is  presumably  related  to  the 
autocorrelation. 

The  only  link  to  an  explanation  may  be  found  from  results  which  were  obtained 
without  a  frequency  window.  Using  the  full  250  Hz  bandwidth  without  frequency  windowing 
led  to  correct  localizations  where  no  localization  was  found  with  windowing.  This  suggests 
that  any  degradation  due  to  sidelobes  of  the  window  function  is  outweighed  by  the 
improvement  due  to  the  effective  larger  bandwidth.  An  additional  improvement  was 
achieved  by  accounting  for  the  phase  shift  in  the  'half  match',  which  led  to  a  considerable 
enhancement  of  the  peak  to  sidelobe  ratio.  These  results  may  also  suggest  that  the  required 
bandwidth  is  not  a  monotone  increasing  function  of  range  and  the  sound  speed  profile 
probably  has  a  distinct  influence  on  that  relationship. 

For  the  iso-speed  environment,  good  localization  is  found  over  the  whole  range  from 
900-7500  m.  At  larger  ranges  the  performance  is  better  than  at  short  ranges  due  to  the  higher 
sidelobes  at  short  ranges.  Footprints  are  slightly  smaller  than  for  the  baseline  results.  These 
results  were  not  found  to  be  very  sensitive  to  the  time  sampling  rate,  in  contrast  to  the 
baseline  results. 
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Figure  9   Autocorrelation  matching  using  a  PE  generated  source 
signal  and  ray  code  templates. (a)  Negative  gradient 
environment  with  a  source  located  at  5500  m  range  and  59.8  m 
depth  and  a  receiver  at  4  0.2  m  depth.  Full  250  Hz  bandwidth 
has  been  used  without  window  and  the  phase  difference  of  the 
autocorrelations  has  been  accounted  for.(b)  Iso-speed 
environment  with  source  at  5000  m  range  and  59.8  m  depth  and 
receiver  at  40.2  m  depth. 
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L  SIFT  MATCHING  RESULTS  (LOW  FREQUENCY) 

The  SIFT  matching  algorithm  that  was  developed  by  Miller  and  others  (Miller  et  al, 
1995),  has  been  modified  to  account  for  bandpass  signals,  to  bring  it  closer  to  a  true 
autocorrelation  matching,  to  improve  the  notching  for  short  time  lags,  to  adapt  it  for  1-  to 
n-point  matching  instead  of  the  fixed  3-point  matching,  to  manipulate  the  amplitude  of  the 
cross  correlation  of  the  source  function  and  the  Hilbert  transform,  and  to  include  the  ray 
amplitude  in  the  matching.  The  original  SIFT  algorithm  is  based  on  a  source  signal  that 
behaves  like  a  low  frequency  pulse  signal.  This  type  of  signal  has  been  proven  adequate  for 
modeling  of  SUS  charges.  However,  it  may  be  less  suitable  to  describe  other  signals, 
especially  signals  which  have  been  high-pass  filtered.  For  this  kind  of  signal  an  n-point 
match,  where  n  increases  with  decreasing  bandwidth,  may  be  more  appropriate  although  the 
matching  may  become  too  complicated  for  large  n.  In  certain  favorable  cases,  tuning  of  the 
6t  factor  in  Eq.  81  may  even  allow  the  use  of  a  3-point  match  when  using  a  signal  that  does 
not  resemble  the  assumed  signal  model. 

The  original  SIFT  algorithm  will  never  give  a  perfect  match  even  if  the  received 
signals  can  be  modeled  as  spikes.  This  is  because  the  cross  correlations  are  stacked  as  an 
array  rather  than  summing  the  cross  correlation  of  individual  arrivals, .  However,  this  result 
is  matched  with  the  sampled  autocorrelation  of  the  received  signal  in  which  the  individual 
arrivals  are  inherently  added  coherently.  As  this  addition  slows  down  the  algorithm,  it  has 
not  always  been  used. 


67 


The  notching  previously  was  based  on  the  arrival  time  difference.  It  did  not  take  into 
account  the  approximated  cross  correlation  of  two  arrivals  containing  two  points  that  were 
shifted  by  6t.  In  the  algorithms  used  here  this  slight  modification,  which  may  be  important 
at  short  lags,  was  taken  into  account. 

Algorithms  for  the  3-point  matching  have  been  implemented.  The  n-point  matching 
algorithm  has  been  considered  but  given  the  time  scale  it  was  not  implemented.  As  the 
source  function  of  the  received  signal  was  known,  the  amplitude  of  the  3-point 
approximation  of  the  cross  correlation  has  been  tuned  to  optimize  the  performance.  For  an 
n-point  SIFT  algorithm,  this  would  be  an  immense  task  and  not  conducive  to  working  with 
real  data. 

The  amplitude  of  the  rays  has  been  brought  back  into  the  SIFT  algorithms  in  contrast 
to  previous  versions.  In  general,  this  improved  the  results  but  only  slightly.  A  few  cases  even 
showed  a  degradation  of  the  results  when  ray  amplitudes  were  included  in  the  algorithm. 

Aside  from  the  parameters  mentioned  above,  the  following  parameters  were  taken 
into  consideration  to  optimize  the  results:  the  offset  6t,  the  sampling  rate,  and  the  bandwidth. 

To  develop  some  intuition  about  the  general  performance,  the  low  frequency  version 
of  the  algorithm  has  been  applied  to  the  negative  sound  speed  and  iso-speed  environments. 
The  results  imply  that  the  performance  of  the  algorithm  for  a  certain  environment,  source 
function,  bandwidth,  and  bottom  depth  is  range  limited.  An  empirical  analysis  suggests  that 
the  area  where  the  algorithm  performs  succesfully  is  below  the  curve 
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d>k    ™  (92) 

\   B 


where  d  is  the  water  depth,  c  is  an  average  sound  speed,  B  is  the  true  bandwidth,  r  is  the 
range  and  k  is  a  constant  that  depends  mainly  on  the  environment.  The  constant  k  ranges 
from  approximately  0.35  to  approximately  0.5  for  the  iso-speed  and  negative  gradient 
environments.  The  approximate  curves  for  the  iso-speed  and  negative  gradient  environments 
are  depicted  in  Fig.  10.  There  is  strong  evidence  that  the  curves  describe  the  areas  where 
the  algorithm  is  likely  to  perform,  but  not  enough  analysis  has  been  done  to  exactly  localize 
the  curves.  The  formula  is  based  on  a  first  order  approximation  of  the  arrival  time 
differences  for  r»d. 

Results  suggest  that  there  may  be  different  regimes  which  govern  the  matching.  First, 
there  are  the  high  bandwidth  cases  which  are  well  below  the  curve  described  by  the 
parameters.  It  is  suggested  that  the  match  consists  of  different  contributions,  most  probably 
energy  and  shape  of  the  arrival  structure,  which  are  influenced  by  the  number  of  notched  out 
values.  For  high  bandwidths,  the  shape  of  the  arrival  structure  prevails.  When  the  bandwidth 
is  reduced,  the  shape  of  the  match  stretches  along  the  lines  of  high  correlation/energy  (most 
probably  corresponding  to  rays  with  dominating  amplitudes).  Reducing  the  bandwidth 
further  causes  the  match  to  split  up  into  several  'matching  points'  which  migrate  along  the 
lines  of  high  correlation/energy.  The  improvements  mentioned  earlier  give  some  relief  but 
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only  marginally  increase  the  performance.  Plots  of  the  number  of  sampled  cross  correlations 
with  lag  values  within  one  pulse  width  do  not  show  any  abrupt  changes  with  decrease  in 
bandwidth. 
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Figure  10   Estimated  performance  limitations  for  the  LF  SIFT  algoritm  for  (a)  the 
negative  gradient  environment  and  (b)  the  iso-speed  environment. 


The  autocorrelation  of  the  received  signal  however  gives  a  clue  to  the  process. 
During  the  reduction  of  the  bandwidth,  a  low  frequency  trend  starts  to  dominate  the  arrival 
structure  of  the  received  signal.  The  individual  arrivals  seem  to  ride  on  top  of  this  low 
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frequency  trend.  In  the  autocorrelation,  the  sequence  of  the  positive/negative  going  parts  of 
the  arrival  structure  becomes  dominant.  This  makes  it  necessary  to  notch  out  about  twice  the 
number  of  points  notched  out  at  higher  bandwidths.  At  these  low  bandwidths,  performance 
is  not  impressive  and  the  footprints  are  large  and  only  approximate  localization  is  possible, 
but  the  sidelobes  are  excessively  large.  It  is  still  unclear  whether  we  could  exploit  this  last 
regime,  especially  in  more  complicated  or  more  noisy  environments.  The  change  in  the 
autocorrelation  of  the  received  signal  indicates  the  limit  up  to  where  the  algorithm  can  be 
pushed  in  the  ordinary  form.  Results  for  both  environments  considered  differ  but  the  trend 
agrees  with  the  theory.  As  the  significant  arrivals  for  the  negative  gradient  environment  are 
more  confined  in  time  than  for  the  iso-speed  environment,  the  limits  for  a  certain  algorithm 
in  this  environment  are  more  severe.  No  low  frequency  synthetic  PE  source  signals  have 
been  used  because  results  of  the  ray  model  and  UMPE  model  are  largely  different  below  10 
Hz  due  to  the  breakdown  of  ray  theory  at  these  frequencies. 

The  size  of  the  footprint  is  found  to  be  dependant  on  range  and  bandwidth.  In 
general,  a  larger  bandwidth  will  give  a  smaller  footprint.  The  numerical  analysis  presented 
here  indicated  that  a  bandwidth  of  250  Hz  at  a  range  of  5  and  7.5  km  was  not  adequate  for 
the  algorithm  to  succesfully  localize  for  the  negative  sound  speed  gradient  environment. 
Consequently,  some  of  the  goals  of  the  project  based  on  these  parameters  had  to  be  dropped. 

Results  show  only  moderate  sensitivity  to  the  offset  parameter  6t,  and  in  some  cases 
a  forced  1 -point  match  (neglecting  the  cross  correlation  of  the  source  function  and  its  Hilbert 
transform)  performs  just  as  well  as  a  3-point  match.  The  optimum  offset  parameter  showed 
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to  agree  well  with  theory.  Using  the  calculated  amplitudes  in  SIFT  instead  of  setting  them 
to  unity  generally  improves  the  results  except  when  the  amplitudes  are  totally  in  error.  Thus, 
with  amplitude  variability  being  a  major  issue  in  real  data,  this  approximation  may  be 
validated,  in  general. 

To  be  able  to  explain  the  good  performance  of  the  algorithm  in  the  Barents  Sea 
experiments,  we  note  that  the  location  of  the  experiment  was  more  than  250  m  deep,  the 
signals  were  essentially  low  frequency,  and  the  low  frequency  variant  of  the  algorithms  was 
used.  The  bottom  was  modeled  as  a  hard  bottom  and  attenuation  was  modeled  as  rough 
surface  scattering,  a  reasonable  approximation  for  this  area.  Both  synthetic  and  real  (SUS) 
sources  were  used.  The  source  range  was  approximately  4.5  km,  the  source  depth  was  15  m 
and  the  receiver  depth  was  170  m  for  most  of  the  analysis.  Applying  the  results  found  in  this 
thesis  suggests  that  the  larger  depth  was  a  crucial  point  in  the  performance  of  the  algorithm 
in  the  Barents  Sea  analysis.  In  addition,  shallow  source  depth  and  large  grid  spacing  in  range 
cause  any  change  in  the  shape  of  the  matching  point  to  be  resolved  in  the  same  resolution 
cell. 

J.         SIFT  MATCHING  RESULTS  (HIGH  FREQUENCY) 

Both  the  1 -point  and  3-point  high  frequency  SIFT  matching  algorithms  were  found 
not  to  be  able  to  localize  in  the  negative  gradient  environment  beyond  4500  m.  In  addition, 
the  3-point  algorithm  was  not  be  able  to  localize  at  ranges  shorter  than  2000  m.  The  1 -point 
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algorithm  localizes  correctly  at  ranges  shorter  than  2000  m  but  shows  high  sidelobes. 
Beyond  4500  m  both  algorithms  show  an  area  of  high  sidelobes  around  the  actual  source 
location.  This  is  only  seen  for  very  high  bandwidths  (>250  Hz  true  bandwidth). 

The  time  domain  arrival  structure  of  the  received  signal  for  the  3-point  match  shows 
a  large  onset  corresponding  to  the  Hilbert  transform  component  of  the  main  arrival.  This 
results  in  an  offset  in  the  autocorrelation  which  may  partly  explain  the  difficulties  with  the 
algorithm.  Using  ray  code  generated  source  signals  shows  good  performance  for  both 
algorithms  even  at  250  Hz  bandwidth.  This  suggests  that  it  might  even  be  a  programming 
error  since  these  results  should  experience  the  same  difficulties. 

Both  the  1 -point  and  the  3-point  algorithms  show  good  performance  for  the  iso-speed 
environment.  Manipulating  the  individual  amplitudes  in  the  3-point  match  does  improve  the 
results  but  only  slightly.  Other  results,  however,  suggest  that  the  add  and  compare  approach, 
discussed  earlier,  may  improve  the  1 -point  and  3-point  matching  results  at  larger  ranges. 
Some  results  suggest  that  both  high  frequency  algorithms  require  even  more  bandwidth  than 
the  low  frequency  version  of  SIFT.  Unfortunately,  this  might  not  be  possible  in  practice.  It 
is  expected  that  there  might  be  some  better  performance  for  the  3-point  matching  algorithm 
in  the  2-4.5  km  range.  The  more  general  n-point  SIFT  algorithm  has  not  been  used  as  it  was 
expected  to  suffer  even  more  from  the  limited  bandwidth,  and  it  was  considered  to  be  very 
difficult  to  tune  the  individual  dt's  and  amplitudes  of  the  components  correctly. 
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Figure  11   Low  frequency  SIFT  results  for  the  negative  sound 
speed  gradient  environment  using  the  3-point  SIFT  algorithm 
and  250  Hz  bandwidth.  The  scattered  match  results  from  working 
beyond  the  performance  limitations,  (a)  Without  add  and 
compare  approach  applied,  (b)  With  add  and  compare  approach 
applied. 
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Figure  12   HF  frequency  SIFT  localization  results  using  a  PE 
generated  source,  located  at  5500  m  range  and  59.8  m  depth  and 
a  receiver  located  at  40.2  m  depth,  for  a  negative  soundspeed 
gradient  environment,  (a)  1-point  SIFT  algorithm,  bandwidth 
50  0  Hz,  no  window  applied.  Add  and  compare  approach  enhances 
the  result,  (b)  3-point  SIFT  algorithm,  500  Hz  bandwidth,  with 
Blackman  window.  Add  and  compare  approach  enhances  the  result. 
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V.    SUGGESTIONS  FOR  FUTURE  RESEARCH  AND  PRELIMINARY 
ANALYSIS  OF  OTHER  MATCHING  ALGORITHMS 


Directions  for  future  research  for  single  hydrophone  Matched  Field  Processing  is  one 
of  the  goals  of  this  research.  Some  of  the  suggestions  for  future  research  may  be  slight 
modifications  and  improvements  to  existing  algorithms.  Those  suggestions  and  their 
motivation  have  already  been  mentioned  with  the  description  of  the  algorithms  in  Chapter 
3  and  the  analysis  and  results  in  Chapter  4.  In  this  chapter,  suggestions  for  more  fundamental 
changes  to  existing  algorithms  and  new  algorithms  are  given. 

A.        SOURCE  DEPENDENCE 

The  dependence  on  the  source  spectrum  has  been  neglected  throughout  this  research. 
If  the  source  is  a  single,  broadband,  pulse-like  transient,  the  specific  details  of  the  spectrum 
may  not  be  critically  important  A  standard  practice  may  be  to  filter  out  some  portion  of  the 
signal  spectrum  over  which  all  the  processing  will  be  done.  An  algorithm  which  is  fairly 
insensitive  to  the  spectral  details  is  then  needed. 

However,  if  the  source  is  composed  of  two  or  more  closely  spaced  pulse-like  events, 
the  separation  of  the  various  multipaths  at  the  receiver  location  becomes  a  formidable  task. 
At  issue  then  is  whether  a  localization  algorithm  can  be  effective  when  utilizing  only 
snippets  of  the  arrival  (i.e.,  portions  of  the  signal  believed  to  contain  multipaths  from  a 
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single  pulse-like  event).  Presumably,  this  is  true  if  enough  single  event  multipath 
information  is  available  for  analysis.  The  degradation  of  the  localization  due  to  multipath 
information  needs  to  be  investigated  further. 

B.         MATCHING  THE  TIME-DOMAIN  ARRIVAL  STRUCTURES 

Earlier  analysis  which  directly  compared  arrival  structures  indicated  that  this 
coherent  approach  is  too  sensitive  to  phase  errors.  The  autocorrelation  matching  technique 
also  suffers  from  phase  errors  in  the  phase  differences  between  separate  arrivals.  This 
suggests  that  we  first  remove  all  phase  information  by  taking  the  absolute  value  of  the 
basebanded,  time-domain  arrival  structure.  In  order  to  emphasize  the  smaller  amplitude  of 
later  arrivals,  it  might  also  be  beneficial  to  then  apply  a  high-pass  filter  in  the  frequency 
domain  to  remove  the  mean  and  low-order  trends  from  the  time-domain,  absolute  value 
arrival  structure.  On  this  signal  you  can  apply  the  standard  autocorrelation  matching  or  just 
apply  a  normalized  cross-correlation  with  a  replica  to  perform  the  matching.  Careful  design 
of  the  high-pass  filter  is  required,  however,  due  to  the  transient  response  of  the  filter  in  the 
time-domain  (The  large  peak  in  the  arrival  structure  due  to  the  initial  arrivals  triggers  a 
large  transient  response).  Figs.  13  (a)  and  13  (b)  show  the  absolute  value  of  the  arrival 
structure  and  the  demoduted  signal. 
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Figure  13   Arrival  structure.   (a)  Absolute  value.   (b)  After 
high-pass  filtering. 
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A  limited  initial  analysis  has  been  done  with  some  success.  As  the  problem  requires 
a  specific  filter  which  creates  a  minimum  of  phase  and  amplitude  distortion  a  low 
coefficient  Butterworth  zero  phase  filter  algorithm  has  been  used.  For  the  complete  range 
in  the  iso-speed  and  negative  gradient  environment,  a  filter  with  2  coefficients  (effectively 
4  coefficients  due  to  the  zero  phase  filter  algorithm)  and  a  cutoff  of  0.3  times  the  Nyquist 
frequency  using  a  sampling  frequency  of  250  Hz  gave  reasonable  results.  (Nb.  The  cutoff 
is  related  to  the  bandwidth  and  the  sampling  frequency).  As  the  parameters  of  the  matching 
algorithm  have  not  been  tuned  to  optimum  performance,  results  may  change  a  little  after 
tuning  of  the  parameters  (e.g.  number  of  notched  lag  values). 

For  autocorrelation  matching,  the  baseline  result  using  a  negative  gradient 
environment  gives  a  footprint  at  5500  m  which  is  80  cm  in  depth  and  approximately  40  m 
in  range.  Using  a  coherent  match,  the  footprint  becomes  slightly  larger  in  range  and 
considerably  larger  in  depth,  approximately  120  cm  in  depth  and  40  m  range.  Sidelobe 
levels  are,  however,  considerably  higher  than  for  standard  autocorrelation  matching  or 
coherent  matching.  Preliminary  results  using  a  ray  code  based  template  show  that  the  direct 
match  (normalized  cross  correlation)  and  the  autocorrelation  may  perform  very  well  for 
some  source  locations,  while  for  other  source  locations  performance  may  degrade  severely 
due  to  high  sidelobes.  All  footprints  were  found  to  be  very  small  in  depth. 
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Figure  14  Autocorrelation  matching  with  high-pass  filtered 
arrival  structures  for  the  negative  sound  speed  gradient 
environment,  (a)  Baseline  result,  for  a  source  at  5500  m  range 
and  59.8m  depth  and  a  receiver  at  40.2  m  depth.  Note  the  high 
sidelobes  and  the  large  footprint,  (b)  Results  using  a  PE 
generated  source  signal  and  ray  code  generated  templates. 
There  are  small  localization  errors  in  range  and  depth. 
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It  has  been  shown  that  the  sample  rate  has  a  significant  influence  on  the  results.  For 
the  iso-speed  environment,  similar  results  have  been  found.  The  footprints  for  this 
environment  are,  however,  larger.  Combining  this  approach  with  a  SIFT-like  1 -point 
approximation  algorithm  did  not  perform  well.  This  approach  only  relies  on  arrival  times 
and  completely  neglects  the  phase  and  frequency  for  making  the  templates. 

C.        PARAMETRIZATION  OF  THE  TIME-DOMAIN  ARRIVAL  STRUCTURE 

Although  originally  not  incorporated  in  the  project,  some  analysis  has  been  done  with 
parametrization  of  the  arrival  structures  and  the  replicas.  As  the  arrival  structures  decay 
rapidly,  it  was  suggested  that  they  could  be  modeled  as  an  auto  regressive  process.  For  the 
analysis,  linear  prediction  coefficients  were  used  and  the  Euclidian  distance  was  used  as  a 
matching  criterion.  Results  of  the  matching  suggested  that  matching  is  possible  without 
having  to  align  the  signals.  Taking  it  a  bit  further  would  suggest  that  the  same  approach 
could  be  applied  to  one  side  of  the  autocorrelation  function. 

A  limited  preliminary  analysis  of  the  results  has  been  performed.  The  scale  used  to 
present  the  results  is  log(l/D)  where  D  is  the  euclidian  distance.  The  matching  works  like 
a  signature.  The  baseline  footprint  is  very  small,  1  resolution  cell  (5  m  x  40  cm).  A  limited 
analysis  of  the  robustness  has  been  done  using  a  random  uniform  20%  disturbance  of  the 
amplitudes  of  the  frequency  bins  of  the  PE  results.  The  match  will  stretch  in  range  to  2-5 
resolution  cells  and  2  resolution  cells  in  depth.  The  algorithm  works  better  for  a  low  number 
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of  coefficients.  At  a  range  of  5500  m  for  the  negative  gradient  environment,  performance 
is  severely  degraded  using  more  than  10  coefficients  which  suggests  that  the  description  of 
the  signals  must  not  be  too  detailed.  No  phase  perturbations  have  been  tested.  It  is  assumed 
that  this  might  more  severely  degrade  the  performance.  Using  a  PE  generated  source  signal 
and  ray  code  generated  templates  the  algorithm  still  localizes  but  is  severely  degraded  by 
high  sidelobes.  This  preliminary  analysis  has  only  been  performed  in  the  4500-6500  m 
range. 
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Figure  15   Results  for  parameterized  matching  for  a  negative 
sound  speed  gradient  environment,  where  arrival   structures 
are  represented  by  7  linear  prediction  coefficients, 
(a)  Baseline  results,  for  source  located  at  5500  m  range  and 
59.8  m  depth,  and  receiver  at  40.2  m  depth.  The  results  is 
hardly  visible  as  the  match  works  like  a  signature  and 
measures  only  one  resolution  cell,  (b)  Results  using  a  PE 
generated  source  and  ray  code  templates.  Note  the  high 
sidelobes . 
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D.        AUTOCORRELATION  MATCHING  OF  THE  FREQUENCY  DOMAIN 
RESPONSE 


The  received  signal  shows  a  particular  signature  in  the  time  domain  as  well  as  in  the 
frequency  domain  due  to  the  multipath  interference.  Note  that  the  absolute  value  of  the  time- 
domain  arrival  structure  described  above  is  related  to  the  autocorrelation  of  the  frequency 
domain  response  through  a  Fourier  transform  by 


GJf)^R{xW(f+X)dl 

— oo 

(93) 


fR(t)R\t) 


e  i2^dt. 


A  frequency-domain  autocorrelation  matching  algorithm  similar  to  the  time-domain 
approach  would  then  produce  an  ambiguity  surface.  The  autocorrelation  matching  can  now 
be  formulated  as 


A=fG^(f)Gpp(f)df  (94) 
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and  we  can  again  notch  out  values  near  zero  "lag"  to  remove  the  mean  and  low-order  trends 
of  the  time  domain.  Note  that  this  technique  is  essentially  equivalent  to  matching  the 
absolute  value  of  the  time-domain  arrival  structures. 

Formulating  the  matching  of  the  autocorrelation  of  the  frequency  response,  however, 
has  reintroduced  the  problem  of  defining  an  absolute  reference  time.  In  the  definition  of  the 
replica  autocorrelation,  it  is  therefore  necessary  to  introduce  a  free  parameter  to  allow  for 
searching  in  the  time  domain,  i.e., 


Gpp{f,x)=fp{t-z)p  \t-x)e  -Wdt 

-oo 

OO 

^Pit')p\t')e-i2^'^dt'  (95) 

—  oo 

oo 

=ea«*fP(X)P\x+JW- 


The  autocorrelation  matching  would  now  look  like 


^)=fGm(f)GPp(fsW-  (96) 


The  values  used  for  the  ambiguity  surface  are  now  the  values  of  A(r)  for  which  A  is 
maximum.  Furthermore,  this  algorithm  applies  to  both  real  and  complex  signals,  an 
improvement  over  the  original  SIFT  design. 
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It  should  be  noted  that,  while  we  developed  this  idea  independently,  it  was  recently 
brought  to  our  attention  that  other  investigators  have  previously  used  an  extremely  similar 
approach  with  reasonable  success  (Nghiem-Phu  et  al.,  1992).  We  believe  they  were  the  first 
to  suggest  using  this  approach  for  this  problem.  However,  they  did  not  apply  any  of  the 
filtering  techniques  we  have  suggested  and  improvements  to  their  original  concepts  may  still 
be  forthcoming. 
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VI.   CONCLUSIONS 

Results  of  autocorrelation  matching  and  approximate  autocorrelation  matching  at  the 
primary  range  of  interest  (5-10  km)  proved  to  be  highly  inconsistent  and  sensitive  to  signal 
mismatch  between  ray-based  methods  and  full-wave  techniques.  Due  to  our  inability  to 
generate  data  in  the  2-4.5  km  range  for  the  negative  sound  speed  gradient  environment,  only 
limited  conclusions  can  be  drawn  from  the  results.  The  autocorrelation  matching  results 
appear  to  be  bandwidth  limited  for  longer  ranges  and  limited  by  the  sidelobe  levels  for 
shorter  ranges.  Results  suggest  that  using  the  full  250  Hz  bandwidth  without  windowing  can 
lead  to  succesfull  localization  at  ranges  of  more  than  5  km.  Footprint  sizes  are  very  small 
in  depth  (<  1  m)  and  moderately  small  in  range  (<15  m).  Taking  into  account  phase 
differences  in  the  autocorrelations  improves  the  peak  to  sidelobe  levels.  Notching  out  values 
near  zero  lag  improves  the  dynamic  range  and  reduces  the  chance  of  invalid  localizations 
(due  to  energy  matching  instead  of  shape  matching)  but  reduces  the  footprint  size.  Results 
for  the  iso-speed  environment  are  much  better,  as  may  be  expected,  and  footprints  are  larger 
than  for  the  negative  gradient  environment. 

The  SIFT  algorithm  is  very  sensitive  and  more  bandwidth  limited  than  the 
autocorrelation  matching.  For  the  low  frequency  version  of  the  algorithm,  a  first  order 
description  of  performance  limitations  is  suggested.  The  high  frequency  versions  of  the 
algorithm  did  not  give  consistent  results  for  the  negative  gradient  environment  in  either  the 
750-1000  Hz  or  the  500-1000  Hz  band  but  proved  to  work  well  in  the  iso-speed 
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environment.  The  high  frequency  versions  seem  to  be  even  more  bandwitdh  limited  than  the 
low  frequency  versions  and  thereby  severly  reduces  the  range  of  applicability.  Results 
suggest  that  localization  can  be  performed  at  ranges  less  than  5  km  for  a  bandwidth  of  250 
Hz  and  a  bottom  depth  of  100  m  and  moderate  sound  speed  gradient  environments.  Using 
a  large  bandwidth  'localization  by  sidelobes'  could  be  performed  in  some  cases  but  with 
minimal  resolution.  Finally,  additional  improvements  may  be  achieved  for  the  negative 
gradient  environment  by  taking  into  account  the  phase  shift  in  the  "half  matching 
autocorrelation  results,  using  no  window  (thereby  using  the  full  bandwidth)  and  extending 
the  range  to  include  the  2^.5  km  bracket.  The  same  approach  may  improve  the  SIFT  results 
where  the  multi-point  approach  may  also  be  considered. 
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APPENDIX  :  SPATIAL  BEAMFORMING  AS  A  FORM  OF  MFP 

The  autocorrelation  matching  in  the  frequency  domain  can  be  related  to  a  generalized 
and  normalized  Bartlett  beamformer.  The  main  difference  is  the  neglect  of  the  off  diagonal 
terms  in  the  correlation  matrix  of  the  frequency  domain  received  signal.  This  is  closely 
related  to  the  fact  mentioned  by  Miller  et  al.  (1996)  that  the  correlation-based  algorithms  are 
likely  to  work  better  on  noise-like  signals  than  on  sinusoidally  signals. 

The  approach  presented  here  is  an  analogue  of  the  spatial  techniques  known  in 
beamforming  and  frequency  bank  techniques  known  in  spectral  estimation.  In  most 
applications  of  Matched  Field  Processing  (MFP),  an  array  of  receiver  hydrophones  is  used. 
The  output  of  the  beamformer  for  the  standard  Bartlett  beamformer  is  then  defined  as 

BB=E{\b(tf}  =weHRwQ=w  HE{p{t)p  H(t)}w0  (97) 

where  Rp  is  the  matrix  of  the  data  p(t),  and  w  is  the  weighting  vector  for  a  particular  look 
direction.  The  data  p(t)  can  be  defined  as 

A(O=«,(0/>/(O+",(O       *=1, #,  (98) 

where  N  is  the  number  of  elements  in  the  array,  n(t)  is  the  noise,  at  (0)  is  the  array  response 
for  a  source  in  the  direction  01  and  s,(t)  is  the  source  function  for  a  source  in  direction  0, . 
Usually  this  approach  is  based  on  the  plane  wave  approximation,  and  the  beamformer  can 
be  refined  to  multi-variance,  multiple  constraint  or  more  exotic  schemes.  The  step  from 
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plane  wave  beamforming  to  MFP  is  made  by  including  the  knowledge  of  the  environment 
and  the  propagation,  usually  by  means  of  a  propagation  model,  in  the  weighting  vectors. 
This  enables  you  to  make  estimates  in  more  dimensions,  usually  a  2-D  range/depth  plane. 
Necessary  for  this  generalized  beamforming  to  work  is  that  the  description  of  the 
environment  and  the  propagation  is  accurate  enough,  the  environment  has  enough  variation 
over  the  distances  of  interest,  and  the  weight  vectors  at  different  points  are  different  enough 
to  form  a  consistent  match.  In  MFP  the  weighting  vectors  are  usually  called  replicas  and  the 
generalized  beamformers  are  called  processors.  The  most  common  are  the  Bartlett,  the 
multi-variance  and  the  multiple  constraint  processors, 

BB=wHR^ 


BMJwHRp'w 


^=c^VHRpxW 


'MCM 


(99) 

1 


where  Wis  the  matrix  of  the  replica  vectors  defining  the  multiple  constraints  around  the  look 
direction.  The  elements  of  c  are  defined  by 


=wH\ 

m         l       m 


(100) 


where  w,  is  the  replica  vector  in  the  look  direction  and  wm  are  the  vectors  defining  the 
constraints. 

The  above  spatial  approach  can  be  applied  to  a  broadband  time  domain  signal  using 
one  hydrophone  instead  of  the  spatial  array.  Two  different  approaches  will  be  shown,  one 
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in  the  time  domain  and  one  in  the  frequency  domain.  We  define  the  preenvelope  (or 
complex  envelope  of  the  signal)  of  the  received  pressure  as 

p(tt)=as(tt)+n(tt)       i=l, JV,  (10i) 

where  a  is  the  complex  response  of  a  hydrophone  to  a  source  at  location  (rs ,  zs )  to  a 
complex  source  at  time  t^.  The  Bartlett  processor  can  be  formulated  as 

Bb=whRml  (102) 

where  w  are  the  replicas.  The  problem  with  this  time-domain  approach  is  the  lack  of  a  time 
reference  for  the  replicas.  With  an  array,  all  the  hydrophone  locations  are  assumed  known 
and  provide  a  reference  in  space.  The  reference  in  time  is  not  so  easy  to  establish.  Unless 
you  can  overcome  this  problem,  this  approach  does  not  seem  viable 

On  the  other  hand,  taking  the  second  frequency-domain  approach  a  proper  absolute 
time  reference  can  be  established.  Define  the  frequency  domain  received  pressure  as 

pif^asif^+rKf).  (103) 

Treating  it  the  same  way  as  the  spatial  array  Rp  can  then  be  defined  by 

R^W)  £*(/)} 

(104) 

=E{d  mH<MH+£(f>£H(f)}- 
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The  generalized  'beamformer'  or  matching  algorithms  can  be  defined  as 

BB=wHRSv.  (105) 

This  approach  does  not  encounter  the  reference  problems  an  approach  in  the  time  domain 
would  have.  This  approach  could  also  be  extended  to  other  types  of  processors.  The  Bartlett 
processor  is  closely  related  to  autocorrelation  matching.  If  we  define  the  autocorrelation 
matching  in  the  frequency  domain  as 


(106) 


{L  \Gju0L  \GnV)? 


it  is  clear  that  this  looks  like  a  normalized  Bartlett  processor  where  the  off  diagonal  terms 
have  been  neglected. 
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